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INTRODUCTION 


The  distortions  of  a  laser  beam  produced  by  density  and  thermal  varia¬ 
tions  In  a  fluid  medium  have  been  the  subject  of  many  Investigations.  In 
general,  such  studies  on  the  effects  of  heat  deposition  from  the  laser  beam 
upon  the  propagation  of  that  beam  can  be  classified  Into  two  groups,  depending 
upon  whether  time  dependence  Is  considered.  Host  available  experiments  are 
conveniently  understood  by  reference  to  theoretical  studies  of  the  gross 

effects  of  thermal  deposition  and  fluid  motion,  which  assume  that  a  steady 

■> 

state  will  be  achieved  for  the  deflection  and  distortion  of  the  laser  beam. 

On  the  other  hand,  several  years  ago  It  was  shown  that  one  of  the  solutions 
for  beam  propagation  was  unstable,  In  the  sense  that  small  distortions  In  the 
beam  would  be  amplified  If  the  power  level  was  greater  than  a  threshold  determined 
by  the  thermal  conductivity  of  the  fluid.  For  the  most  part,  those  Investiga¬ 
tions  that  do  assume  some  time  dependence  for  the  beam  pattern  were  concentrated 
on  studies  of  the  growth  of  small  amplitude  perturbations  from  an  Initially 
uniform  beam. 

The  work  described  In  this  report  was  motivated  by  the  discovery  of 
Instabilities  In  the  system  of  equations  describing  electromagnetic  wave 
propagation  and  fluid  dynamics.  Initially,  It  appeared  desirable  to  attempt 
an  exploration  of  the  extended  development  of  the  Instabilities  noted  for  a 
uniform  beam.  The  basic  equations  are  discussed  In  Part  I,  and  In  Part  II  the 
linearized  stability  analysis  Is  presented  along  with  an  evaluation  of  the 
threshold  for  these  growing  waves.  To  follow  the  growth  of  the  disturbances, 
computer  studies  were  undertaken.  In  the  course  of  these  studies  It  became 
apparent  that  there  was  some  merit  to  Introducing  a  new  concept  to  judge  the 
value  of  an  algorithm  for  computing  the  solutions  of  a  system  of  partial 
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differential  equations.  This  concept  was  called  "utility11*  and  Is  discussed 
In  Part  III*  along  with  several  examples  of  Its  application  to  simpler  partial 
differential  equations.  The  advantage  of  this  concept  Is  that  It  Is 
relatively  easy  to  apply  to  complicated  systems  of  partial  differential 
equations,  whereas  the  stability  concept  leads  to  a  very  complicated  procedure 
for  deciding  on  the  value  of  a  numerical  routine.  In  Part  IV  are  presented 
the  results  of  a  calculation  of  beam  distortion  for  a  very  high  Intensity 
pulse  propagating  through  air  for  several  kilometers.  Analytical  arguments 
are  advanced  which  suggest  that  the  qualitative  features  of  the  distortions 
are  correct,  which  lends  credence  to  the  computer  output. 

Speed  and  memory  size  In  a  computer  place  certain  restrictions  on  one's 
ability  to  Investigate  phenomena  In  the  laser  beam  problem.  In  the  attempt  to 
calculate  distortions  of  the  type  predicted  by  the  linearized  Instability 
analysis*  cylindrical  symmetry  was  Imposed  on  the  problem  In  order  to 
facilitate  the  computer  calculation.  Had  this  not  been  necessary,  or  had  some 
other  Independent  variable  been  eliminated  rather  than  the  angle  about  the 
beam  axis,  much  more  pronounced  evidence  of  beam  and  fluid  Instabilities  would 
likely  have  been  observed  for  substantially  lower  powers,  powers  that  may  be 
achievable.  Arguments  supporting  this  proposition  are  contained  In  Part  V. 

Part  VI  of  this  report  contains  an  analytic  discussion  of  beam  bending 
and  thermal  blooming  for  a  slab  beam  propagating  through  a  wind.  A  formula 
Is  derived  which  provides  for  the  transition  between  two  regimes  In  which 
conduction  and  forced  convection,  respectively,  dominate  the  dissipation  of 
heat  deposited  In  the  medium  from  the  laser  beam.  This  formula  appears  to  be 
useful  for  the  analysis  of  several  experiments. 
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BASIC  EQUATIONS 


When  an  Intense  laser  beam  propagates  through  a  fluid,  many  Interesting 
phenomena  take  place.  This  laser-fluid  system  can  be  described  by  a  macro¬ 
scopic  model  which  Involves  Maxwell's  equations,  the  Navler-Stokes  equation, 
an  energy  conservation  equation,  and  the  continuity  equation  for  fluid  motion. 
These  equations  which  describe  the  behavior  of  Intense  electromagnetic  beams 
and  the  associated  sound  and  thermal  fluctuations,  are  coupled  by  stimulated 
Raman  scattering,  electrostrlctlon,  the  high  frequency  Kerr  effect,  absorption 
heating,  and  the  density  and  temperature  dependence  of  the  dielectric  constant. 
In  this  paper,  a  systematic  discussion  Is  presented  for  an  Intense  laser  beam 
propagating  through  air,  which  has  a  negligible  Kerr  constant.  If  the 
frequencies  are  outside  the  Raman  scattering  range,  the  Instabilities  are 
primarily  caused  by  optical -accoustic  coupling  of  the  laser  beam  and  the  gases. 
These  effects  are  of  long  duration  compared  to  those  of  self- focusing.  As  the 
beam  passes  through  air,  the  Intensity  profile  Induces  a  nonuniform  temperature 
gradient  transverse  to  the  propagating  direction  of  the  beam,  due  to  the  energy 
absorption  from  the  beam.  This  thermal  non-equilibrium  and  electrostrlctlon 
together  cause  the  generation  of  a  density  gradient  and  hence  a  sound  wave. 
These  density  changes  react  back  on  the  Incident  beam  through  changes  In  the 
dielectric  constant.  A  detailed  mathematical  description  of  the  model  adopted 
here  will  be  presented. 

Maxwell's  equations  for  the  electric  field  vector  E  In  a  charge  and 
current  free  region  can  be  employed  to  yield  the  wave  equation  for  the 
electric  field 


V  X  (V  X  E) 


(1) 


i  h 
?  it* 


where  the  permittivity,  e,  is  taken  to  be 

c  -  eL  +  t2  ?  (2) 

Here  and  throughout  this  paper,  Heavysl de-Loren tz  units  are  used  for  electrical 
quantities.  These  units  are  sometimes  called  rationalized  gausslan  units.  The 
time  average  In  (2)  Is  to  be  taken  over  a  time  large  compared  with  an  optical 
period.  Many  such  time  averaged  terms  will  be  encountered  below  and  In  all 
cases  an  average  over  "several"  optical  periods  Is  Intended.  Such  average 
quantities  may,  of  course,  still  be  functions  of  time,  but  will  vary  slowly  on 
the  time  scale  of  an  optical  period. 

For  an  Isotropic  medium  free  of  charge,  the  electric  displacement 
vector  Is  parallel  to  E  and  Is  divergenceless.  Under  these  conditions  one, 
therefore,  has 

7  •  (eE)  -  0  (3) 

From  (3)  one  obtains 

+  ■+ 

V.  •  E  -  -  E  •  V  (log  e). 

Substituting  this  expression  In  (1)  produces 

V2  E  +  V  (E  •  V  log  e)«  X  -ty)  (4) 

cc  3t 


Since  processes  Involving  acoustical  and  thermal  effects  are  considered 
here,  the  change  In  e  over  an  optical  wavelength  s  small  compared  with  the 
changes  over  a  typical  acoustical  length.  Thus  the  term  Involving  Ve  will  be 
neglected  In  (4),  leaving: 

cVe**^  (eE).  (5) 

ar 
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When  processes  involving  electronic  states  of  the  molecules  of  the  fluid  are 
Important,  the  Ve  term  must  be  retained  in  the  wave  equation  for  E. 

The  permittivity  e  is,  in  general,  a  function  of  the  mass  density  p 
of  the  fluid  and  the  temperature  T.  In  fluids  will  be  small  except  for 
such  large  temperatures  that  there  is  significant  population  of  vibrational 
modes  of  the  molecules.  When  anharmonic  effects  In  the  vibrational  spectrum 
become  important,  careful  attention  must  be  given  to  the  dependence  of  e  on  T 
at  constant  density.  For  a  considerable  range  of  density  and  temperature. 


••  ¥). 


however,  i-£t]  w  0,  so  that  e  can  be  taken  to  be  a  function  of  p  only.  The 
fP 

dependence  of  e  on  density  can  be  approximated  by  the  formulas  resulting  from 
Lorentz-Lorenz  local  field  theory.  This  approach  leads  to  the  formula 

P(l)TB(e'1)iet2)*  f  <«L  -'><«-  «L>  (6> 

The  second  term  on  the  RHS  of  (6)  Is  usually  extremely  small  because  (e  -  eL), 
the  nonlinear  piece, Is  small.  Furthermore,  particularly  in  gases,  (eL  -  1)  is 
extremely  small.  For  the  situations  considered  In  the  present  analysis,  both 
these  factors  are  small,  so  that  (6)  Is  conveniently  simplified  to 


(6a) 


This  Clausius- Mo ssottl  relation  will  be  used  later  in  the  detailed  analysis  of 
the  laser-fluid  equations. 


The  equation  for  the  fluid  motion,  the  Navier-Stokes  equation,  is 


P9j 


(7) 


where  v  Is  the  velocity  of  a  material  element  of  the  fluid,  g  Is  the  gravita¬ 
tional  acceleration  vector,  o^  Is  the  viscosity  stress  tensor, ^  ajj  Is  the 
interaction  stress  tensor  coupling  the  electromagnetic  field  to  the  fluid. ^ 
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The  time  derivative  ^  Is  the  "material"  derivative;  l.e.,  one  that  follows 
the  motion  of  a  material  "particle"  of  the  medium  relative  to  a  fixed  coordinate 
system.  This  derivative  can  be  expressed  In  the  form 

>•  <8> 


The  viscosity  stress  tensor  Is  given,  to  a  first  approximation,  by  the 
linear  terms  of  an  expansion  In  powers  of  the  viscosity  coefficients: 

3V, 


/3Vj  3V,  \ 

*1J  *  +  ^*1/  ^ 


_ k  . 

3\  « 


(9) 


where  n  Is  the  coefficient  of  shear  viscosity  and  n'  Is  the  compresslonal 
viscosity  coefficient.  Sometimes  n  and  n1  are  called  the  "first"  and  "second" 
viscosity  coefficients  and  were  previously  denoted  \i  and  ? ,  respectively.  In 
the  notation  used  by  Stokes.  The  coefficient  n'  Is  occasionally  called  the 
dllatlonal  viscosity  coefficient. 

The  Interaction  stress  tensor  for  the  electromagnetic  field  and  the 
fluid  will  be  taken  to  be 

-2 


'  ll  ’  'P41j  "  T  [E  ’  p(fp)T]  5U  +  e¥j 


(10) 


A  derivation  of  this  tensor,  valid  for  static  electromagnetic  fields  Is  given 
on  page  67  of  Reference  [2].  In  a  vacuum,  this  expression  becomes  the  familiar 
Maxwell  stress  tensor 

°'ij M  EiEi '  ?e2  4ij  • 

This  tensor  Is  not  strictly  correct  for  optical  fields.  Expression  (10)  results 
from  a  derivation  with  an  Isothermal  constraint.  A  similar  derivation  with  an 
adiabatic  constraint,  actually  an  Isentroplc  constraint,  gives  the  same  result, 
except  that  the  partial  derivative  at  constant  temperature  Is  replaced 
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by  ,  the  derivative  at  constant  entropy.  The  difference  In  these  two 
constraints  Is  contained  In  the  thermodynamic  relation 


/9e\  .  „/ae\  .  T/3e\ 
P\3p/T  \3T/ 


p  YCp 


where  Cp  Is  the  specific  heat  per  unit  mass  at  constant  pressure,  y  Is  the 
ratio  of  specific  heat  at  constant  pressure  to  that  at  constant  volume,  @  Is 
IpPMhermal  expansion  coefficient  ■  *  p  »  and  vs  the  isentroP*c 


velocity  of  sound ^ j  • 


The  difference  term  In  (11)  Is  very  small  for  reasons  that  have  already 

_ ,  /3e\  „  TL.  J. _ J  - _ _ _ J  -  — _ L.1  ..  1  Tlui. 


been  discussed 


,  «  0.  The  term  In  square  brackets  Is  roughly  Thus, 


for  the  present  purposes.  It  will  suffice  to  use  the  "Isothermal"  stress  tensor 
given  In  (10).  Actually,  neither  constraint  Is  strictly  valid,  but  corrections 
would  be  small  and  would  necessitate  a  detailed  examination  of  fluid  boundary 
layers  and  the  explicit  mechanisms  of  heat  deposition  In  the  control  volume. 

The  pressure  P  occurring  In  (10)  Is  the  thermodynamic  pressure.  This 
pressure  Is  not  Identical  with  the  mean  pressure,  Pffl  =  -  -j  +  a^],  but  Is, 
rather,  the  thermodynamic  Intensive  variable  that  enters  Into  the  equation  of 


state  of  the  fluid: 


P(p.T) 


Combining  equations  (7),  (9),  and  (10),  one  obtains 


M  -  *  +  %s  +  W 


where  the  electrostrlctlve  force  density  ffi  Is  obtained  by  substituting  (10)  In 
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%s  ■  <  -  7  ’to2)  +  ’[t  o(^)t]  +  C(E  -  V)  E  > 

■  <  -  ?  E2[(|f)T  »P  +  (|f)p  "]  -  e(E  -  V)  f  + 

+  1  e2(|?)t  +  p’[t  (|?)t]  +  p(E  •  ’)  E  > 


The  time  averages  here  are  taken  over  several  optical  periods  and  the  relation 
(3),  V  •  (eE)  ■  0  has  been  used  to  get  (14) y 

The  viscous  force  density  fvjsc.  appearing  In  (13),  Is  obtained  by 
substituting  (9)  In  (7): 

Tyisc  -  n  V2  7  +  (n  +  n’)  v  (v  •  v)  (15) 

where  terms  Involving  Vn  and  Vn'  have  been  dropped.  In  other  words,  n  and  n* 
have  been  taken  to  be  constant  throughout  the  fluid.  This  Is  a  good  approxlma 
tlon  for  gases  where  n,  for  example.  Increases  slightly  with  temperature. 
Models,  based  on  the  Lennard-Jones  potential  energy  function,  suggest  that 
n  -  A  .  In  liquids,  on  the  other  hand,  the  viscosity  decreases  rather 
strongly  with  Increasing  temperature,  so  that  gradients  of  the  viscosity 
coefficients  could  have  some  small  effect. 

Putting  (14)  and  (15)  In  (13),  one  can  write 

P  jft  ’  P9  -  7P  +  P^]  •  7  ^  *  e 

+  n  v2  v  +  (n  V)  v  (v  •  v)  (16) 
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The  equation  of  continuity  of  matter  In  the  fluid  Is 

If  ♦  V  •  <pj)  -  o 

The  general  equation  for  heat  transfer  In  the  fluid  Is 

pT  Dt  m  ♦n  '  v  •  5 


(17  ) 


(18a) 


where  s  Is  the  entropy  per  unit  mass  In  the  fluid,  4>n  Is  the  viscous  dissipation 
function,  and  q  Is  the  total  energy  flux  vector.  The  flux  q  can  be  divided  Into 


several  parts: 


5  "  Votf  +  *rad  *  F  (18b) 

where  F  is  a  "material"  Poyntlng  vector  to  be  discussed  below,  and  qcorKj  and 
qrad  are  the  parts  due  to  conduction  and  thermal  radiation.  In  particular, 

=  -  KV  T  (18c) 


qcond 


s  -  <V  T 


where  the  thermal  conductivity,  <  ,  a  function  of  p  and  T  Is  defined  by  this 
relation.  The  thermal  radiation  flux,  qrad,  can  be  approximated,  for  small 
temperature  differences,  by  Newton’s  law  of  cooling: 

v’  Vad  *  P  cv  *  (T  ’  V  *18d) 

where  C  Is  the  heat  capacity  per  unit  mass  at  constant  volume',  T  -  T  Is  the 

*  °  B 

local  temperature  excess,  and  q  Is  a  radiation  coefficient  Introduced  by  Stokes. 

Because  this  effect  Is  extremely  small  In  the  parameter  regime  of  the  present 

analysis,  this  thermal  radiation  term  will  be  dropped. 

The  viscous  dissipation  function,  $  ,  appearing  In  (18a)  Is  given  by 


E  °1j  9Xj 


(l8e) 


Combining  (18a),  (18b).  (l8c).  and  dropping  the  term  shown  In  (18d),  one  obtains 
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The  LHS  of  (19)  can  be  re-expressed,  and  the  entropy  eliminated,  by 
using  the  thermodynamic  relation 

pT  Dt  “  p  Cv  Dt  "  3  fit  (2°) 

where  the  thermodynamic  parameters  entering  here  have  all  been  defined  above. 
Actually  (20)  Is  obtained  only  If  one  neglects  the  dependence  of  the  state  of 
the  fluid  on  the  electric  field.  Strictly  speaking,  the  entropy  Is  a  function 
of  p,  T,  and  E.  As  shown  on  page  51  of  Reference  T2],  one  has  the  approximate 
expression 


s  «  so(T,p)  + 


(21) 


Since  Is  small,  as  discussed  above,  It  Is  consistent  to  drop  the  very 


small  second  term  In  (21)  and  thus  get  (20). 


Combining  (19)  and  (20),  one  then  gets 


PC, 


DT 

fit 


cv  (Y  •  1) 

3 


ft 


Vv 


(k  V  T)  -  V  •  F 


(22) 


The  term  V  •  F  In  this  equation  Is  associated  with  a  model  for  the  absorption 
of  electromagnetic  energy  by  the  fluid.  A  linear  absorption  coefficient,  a, 
for  the  deposition  of  electromagnetic  energy  In  the  fluid  Is  to  be  Introduced. 
This  coefficient  Is  taken  to  be  Independent  of  the  frequency  of  the  electro¬ 
magnetic  radiation,  so  the  model  will  not  be  valid  for  frequencies  near  the 

resonance  lines  of  the  molecules  In  |he  fluid.  In  this  model,  the  electric 

-  7  2 

field  will  be  damped  by  a  factor  e  where  z  Is  the  direction  of  propagation 

~2 

and  the  energy  deposited  In  the  medium  Is  taken  to  be  ac^c  E  =  aP^  where  P|_ 

o 

Is  the  laser  Intensity  In  (ergs/ sec) /cm1 .  It  Is  this  term  that  the  divergence 
of  F  represents  In  (22).  This  "material"  Poyntlng  vector  F  Is  understood  to 


be  the  time  average 

F«<Ft-?0>  (23) 

where  FT  =  cE  x  B  Is  the  total  Poyntlng  vector,  whereas  FQ  Is  the  Poyntlng 
vector  that  would  result  from  the  same  fields  In  the  absence  of  absorption. 
The  time  average  Is  taken,  as  usual,  over  several  optical  periods. 


In  order  to  be  consistent  In  the  Introduction  of  this  absorption  model, 
certain  terms  should  be  added  to  the  wave  equation  and  to  the  Navler-Stokes 
equation.  Specifically,  the  term  ac  should*  be  added  to  the  RHS  of  (5) 

•f  *4* 

and  a  vector  C  should  be  added  to  the  RHS  of  (16).  The  counter  term  C  Is 
given  by 


(24) 


This  Is  strictly  a  counter  term  whose  purpose  Is  to  remove  the  a  dependence 
from  the  Navler-Stokes  equation.  No  physical  significance  Is  attached  to  this 
term. 


Making  these  modifications  and  collecting  the  operative  equations  for 
convenience,  the  following  set  of  nonlinear  coupled  partial  differential 
equations  Is  obtained  for  description  of  the  macroscopic  representation  of 
the  laser-fluid  system. 


WAVE  EQUATION 


cVE  =  -V(eE) 

ar 


e  *  e. 


+  e„  V 


(25a) 

(25b) 

(25c) 


4* 

Instead  of  adding  this  term  explicitly,  the  same  effect  can  be  obtained  by 
considering  the  permittivity,  e,  to  be  complex  and  frequency  dependent. 

The  present  procedure  Is  used  to  avoid  the  logical  difficulty  of  putting 
frequency  dependence  In  the  space- time  wave  equation. 
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NAVI ER- STOKES  EQUATION 

p  ^  -  p5  -  vp  ♦  ^  v  [e2  p(H)t1  -  i  e2  v  e  + 

+  n  v2  V  +  (n  +  n')  V  (7  •  v)  +  C 
T «  1  0  F 

tT*pWr3c 


HEAT  TRANSFER  EQUATION 


C„(T  -  1) 


•Sr  ft  - «  -  *h  *  T  *  <'"T>t“^E 

3,<  r  F*i  avA  svk  1 8vi 

♦n  H  °1J  aHEJ  '  [y«7  +  SN  j  +  n'  8*k  6ijJ 


FLUID  CONTINUITY  EQUATION 


||  +  V  •  (pv)  «  0 


EQUATION  OF  STATE 


P  *  P(p»T) 


(26a) 

(26b) 

(27a) 

(27b) 

(28) 

(29) 


Due  to  the  complexity  of  this  set  of  equations,  It  Is  not  possible  to  obtain 
exact  solutions  analytically.  Only  the  linearized  solutions  have  been 
discussed. In  the  next  section  a  linearized  analysis  of  this  set  of 
equations  Is  presented  and.  In  subsequent  sections,  the  nonl Inearl  ties  will 
be  considered. 
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PART  II 


LINEARIZED  ANALYSIS 

Linearized  analysis  Is  a  standard  pertubatlon  technique.  In  this 

scheme.  It  Is  assumed  that  each  of  the  dependent  variables  In  the  problem  can 

be  expressed  as  the  sum  of  Its  slowly  varying  zeroth-order  component  and  a 

(6) 

small  first-order  correction.  In  this  way,  a  set  of  linear  equations  for 
small  disturbances  Is  obtained.  This  approach  to  the  analysis  of  the  laser- 
fluid  system  was  first  Investigated  by  Brueckner  and  Jorna.  In  the  present 
approach,  two  variables  are  used  to  describe  the  perturbed  electromagnetic 
field,  one  for  the  component  of  the  field  which  Is  vibrating  In  phase  with  the 
primary  beam  and  one  for  the  component  out  of  phase.  In  this  way,  the  four- 
photon  coupling  Induced  by  periodic  fluctuations  In  the  dielectric  constant 
can  be  Included.  This  coupling  was  not  Included  In  the  original  formulation 
given  by  Brueckner  and  Jorna.  The  dispersion  relation  for  these  linearized 
equations  has  been  evaluated  and  Is  more  complicated  In  structure  than  that 
presented  by  Brueckner  and  Jorna.  For  propagation  through  air,  however,  the 
numerical  differences  are  minor.  The  wave  with  the  largest  growth  rate, 
resulting  from  resonant  Interactions  between  scattered  electromagnetic  waves 
and  the  thermal  wave,  propagates  almost  perpendicularly  to  the  laser  beam. 

The  direction  Is  such  that  the  change  In  frequency  of  the  scattered  electro¬ 
magnetic  wave  and  the  frequency  of  the  thermal  wave,  which  Is  zero,  are 
approximately  the  same. 

-A,4^fl1igd^analys1s  will  now  be  presented.  Separating  each  quantity 
Into  a  zeroth-order  and  a  first-order  perturbation,  one  has 
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(30a) 


E*E<0)  +  E(l) 

P  »  Pq  +  P^ 

T*T0  +  T1 

p“po  +  Pi(l0T  ♦  Ti(lf)p  ■  "o  ♦  \  pi  + 1- 6  po  Ti 

E  *  [£l(o)  ♦  pi(§r)T  +  ti(st^p]  + 

*  [e2(0)  +  pl(a/)T  ♦  Tl(ar)J  (E(0)  +  E(1))2 

"  E0e+  E1 


(30b) 

(30c) 

(30d) 

(30e) 


where  p0,  Tq,  Pg,  e^gj,  and  e2(gj  are  taken  to  be  constants  describing  the 
unperturbed  medium,  v$  Is  the  Isentroplc  velocity  of  sound  In  the  unperturbed 
medium,  y  Is  the  ratio  of  specific  heats,  $  Is  the  coefficient  of  thermal 
expansion  for  the  fluid,  and  the  zeroth  and  first  order  dielectric  constants 
are  given  by 


Eoe  =  El(0)  +  e2(0)  e0 


(30f ) 


[(£). Pl'frX 7]' 


+  2  e 


2(0)  fc(0)  E(l) 


(30B) 


The  indicated  averages  are  taken  over  several  optical  periods  and  the  subscript 
zero  on  the  partial  derivatives  Indicates  that  they  are  to  be  evaluated  at  the 
density  and  temperature  of  the  unperturbed  fluid. 


When  equations  (30)  are  substituted  Into  (25a)  and  the  zeroth-order  and 
first-order  terms  are  separated,  one  gets 


-14- 


1  voe  E(0) 


nbr  •  \  vc 

7E0)-?  H 


A*< 

c  §1 


oe  E(0) 


(31a) 


1  \eoe  E(l))  o  E(1)J  _  1  3  Vel  E(0) 


<’)  «  91  c2  3t2  ' 

There  would  have  been  an  additional  small  term  -$[Eq  •  V  log  e]  on  the  RHS  of 

(31b),  had  equation  (4)  not  been  simplified  by  dropping  Ve  to  get  (14).  Since 

a  3(^E0) 

o  Is  treated  as  a  first  order  quantity,  a  term  — - - tt—  was  also  dropped 

2c/e 

from  the  RHS  of  (31b)  oe 

The  primary  electromagnetic  wave  Is  chosen  to  be  a  linearly  polarized 
plane  wave  propagating  In  the  z-dlrectlon: 


(31b) 


:(0)  ’  I  ey  E0  *  *  +  c-c' 

-  f  z 

■  ey  | Eq |  e  cos  (u^t  -  kLz  +  6) 


(32a) 


where  c.c.  stands  for  "complex  conjugate"  and  6  Is  the  phase  of  the  complex 
constant,  Eq.  With 

-  kLC 


n.  • 

oe 

equation  (32a)  gives  a  solution  of  (31a),  except  for  negligible  terms  of 


(32b) 


order  a  . 

The  first  order  correction  described  by  (27b)  Is  taken  to  be  of  the 


+  +  +  +  a 

-*■  «  f  1  (w+t  -  k.  •  x  +6)  -1(b)  t  -  k  *  x  +  6)1  "  1 

E0)'Iey[fe  +ge  -  •  Je 


+  C.C. 


(33a) 
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where  f  and  g  are  complex  amplitudes,  6  Is  the  phase  of  the  primary  beam,  and 


ta )+  =  +  W 

(33b) 

*  ■+ 
k+  5  k|_  ez  +  k 

(33c) 

k+2  *  kL2  +  k2  +  2  k3  k|_ 

- 

(33d) 

It  Is  presumed  that 

'I 

cd  «  fa)j_ 

(34a) 

and  that 

7T  + 

A 

A 

• 

(34b) 

Under  this  assumption.  It  follows  that  ez  *  k+is  approximately  equal  to  k^.  One 
can  easily  check,  now  that  the  factor 


appearing  in  (33a)  1c  consistent  with  the  term  containing  a  In  (31b)  to  a  very 


good  approximation. 


both  have  the  same  damping  factor 


Because  of  the  mathematical  complexity.  It  Is  convenient  to  express  the 

first  order  electric  field  In  another  form.  Equation  (33a)  Is  rewritten: 

a  z 

E(i )  »  ey  e"^  [Ej  cos^t  -  kj_z  +  6)  +  Ej  sin^t  -  kLz  +  6)]  (35) 


where 


Ej  =(^)  e^wt  -  *  *  *>  +  c.c.  - 

4  4 

■  |f  +  g|  cos(wt  -  k  •  x  +  6' ) 


(35a) 


E'j  =  1  *(f§3-)e1(uJt  '  1  *  *>  ♦  c.c.  = 
■  |f  -  g|  cos(wt  -  k  •  x  +  6") 


(35b) 
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where  6'  Is  the  phase  of  the  complex  number  (f+g)  and  6"  is  the  phase  of 
1(f-g).  Since  <a  «  w^,  the  amplitudes  Ej  and  E^'  are  slowly  varying  functions 
of  time  compared  with  the  rapidly  varying  primary  beam. 

Using  (32a)  and  (35)»  the  time  averages  appearing  In  equations  (30f)  and 
(30g)  are  easily  computed  and  found  to  be: 


■(0) 


“  ?  lEn 


2  -c*z 


(36  ) 


E(0)  E(1)  -  \  | Eq |  E\  e"aZ  -  \  |E0|  |f+g|  e‘“z  cos(ut  -  k 
*  {  |E0|  (f+g)  -  1  ■  *>  e-“z  +  c.c. 

The  first  order  sound  and  thermal  waves  are  taken  to  be 


x  +  6') 


(37  ) 


(?;)■?(?!)  e1(“‘ -  ♦  c.c.  (38) 

where  p1  and  T'  are  complex  amplitudes.  The  frequency  to  and  wave  vector  k 

appearing  in  (38)  are  the  same  as  those  shown  In  (33b)  and  (33c).  These  waves 
+ 

at  frequency  +  (toL  -  u.^)  are  often  called  thermal  Rayleigh  waves  or  Brillouln 
waves. 

The  particle  velocity  v  can  be  taken  to  be  a  first  order  quantity  for 

the  case  where  the  unperturbed  medium  Is  stationary.  For  this  reason  and 

-+•  -+• 

because  the  unperturbed ‘medium  is  taken  to  be  uniform,  so  that  V  pQ  and  V  TQ 
are  zero,  there  Is  no  difference  In  first  order  between  and  the  partial 
derivative  Thus,  all  of  the  material  derivatives  reduce  to  ordinary 

partial  derivatives  and  there  Is  no  distinction  between  Eulerian  and  Lagranglan 
representations  In  the  first  order  calculations.  The  continuity  equation  (28) 
becomes 


-17- 


(39) 


9Pl  -+■ 

ir  +  po7'  v"° 


and  the  Navler-Stokes  equation  (2ba)  becomes 
PoI'Cri  ^Mn  +  n')  v(v-  v)] 


■  \  E(0>  7e  +  * 


V  +  f 

=  -  CVP]  +  0PO  VT-j]  + 

+  [i  p(i)T  7  £2  +  ^i-11  7  £]  * E  (40) 

The  Clauslus-Mossottl  relation  (25c)  and  equation  (30d)  have  been  used  to 
obtain  (40). 


Now  (39)  and  (40)  can  be  combined  In  such  a  way  as  to  remove  v  from  the 
equations.  To  accomplish  this,  (39)  is  differentiated  with  respect  to  time  and 
the  divergence  of  (40)  is  taken.  The  quantity  V  •  v  =  V  •  |^-  Is  then 
eliminated  between  the  resulting  equations  to  produce: 


[V2Pl  +  Bp0  VT-j]  = 
7e+e] 


U  -  DU  +  2)  2 

3 


E(0)  E(l) 


(41) 


where  the  Clausius-Mossottl  relation  has  again  been  used  and  the  extremely  small 

1  o  A 

term,  j  (e  -  l)  aE  ez  •  Ve.has  been  dropped  along  with  other  second  order 
terms.  The  gravity  term  is  also  dropped  In  (41)  because  it  Is  of  no  interest  in 
the  present  considerations.  The  counter  term  £  has  performed  the  task  for  which 
It  was  Introduced.  The  net  effect  of  the  term  is  the  instruction  that  the 
Laplacian  ectlng  on  should  not  act  on  the  damping  factor  e”aZ. 
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For  convenience,  the  following  notation  Is  Introduced: 


2  V 

u  =  .  so  that 

u  =  Isothermal  speed  of  sound 
In  the  unperturbed  medium 


N  =  IgQjLaD 

p0 

o 

In  this  notation,  V  e  becomes,  to  first  order, 

»*•*■»•  [(|)T  *>1  ♦  (lf)fl  "l]  *  +  bA1 

and  equation  (41)  becomes 


(42a) 


(42b) 

(42c) 

(42d) 


(43) 


^l-N^).UVP,t6pA] 


♦  *  "A,)  +  <Se+2)viW(7j'] 

(44) 

To  first  order,  the  thermal  transfer  equation  (27a)  becomes 

p0Cv  at“  "  r~Ws  **l  +  ac’^oe  |_E(0)  +  2E(0)E(1)J  * 

Upon  substitution  of  quatlons  (32a),  (35),  (36),  (  37)  and  (38)  Into 
equations  (31b),  (44),  and  (45),  the  Fourier  transform  of  the  set  of  linearized 
equations  Is  obtained  and  the  corresponding  Fourier  amplitudes  are  related  by  a 
set  of  four  simultaneous  linear  equations.  Two  of  these  equations  result  from 
(31b)  because  two  different  frequencies  enter  Into  that  equation.  The  constant 
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term  ac/e^  E(q)  *s  dr°PPed  from  because  It  corresponds  to  a  uniform 
linear  Increase  In  the  temperature  of  the  medium  which  Is  of  no  Interest  In 
the  dispersion  relation.  The  damping  factor 

-f2 


divides  out  of  (31b),  the  first  order  wave  equation,  but  the  factor  e 


■aZ 


entering  (44)  and  (45)  with  the  terms  does  not  divide  out.  There  are 

two  alternatives:  One  could  go  back  and  try  to  Insert  a  double  damping  factor 


In  (38),  but  this  would  be  rather  artificial.  Instead,  since  a  Is  small,  the 
damping  factor  Is  merely  replaced  by  1  In  (44)  and  (45).  Thus  the  considera¬ 
tions  below  pertain  only  to  values  of  z  small  enough  that  the  decreasing  ampli¬ 
tude  of  the  electromagnetic  waves  does  not  significantly  alter  the  equations 
for  thermal  and  sound  propagation.  The  four  equations  relating  f,  g,  p1,  and 
T'  are: 


free  A-**]  fplW].* 
♦  [i  «o4}'  *  [?  <v+] r  ■ 0 

+  [l  AEo“-2]  o'  +  [j  BEo“-]  T'  ■  0 


(46) 


[ac^oeEoJ  f  +  »‘^eEo]  8  +  [iJtT“  “J  »'  + 


(47) 


Eok2] f  E/] ,  ♦ 


+  p-1N<Dk2-(u2  -  AE2k2j  p'  + 
+  [(-u2SP0  +  BE|)  k2]  t-  .  0 


Throughout  these  equations,  EQ  has  been  written  In  place  of  |EJ.  Since  only 
|t'Q|  appears,  the  phase  6  does  not  enter  these  equations.  Therefore,  without 
loss  of  generality.  EQ  will  be  considered  to  be  real  and  positive.  The  con¬ 
sistency  condition  for  these  four  equations  Is  the  vanishing  of  the  determinant 
of  the  coefficients  of  f,  g,  p'  and  T'.  Since  in  (48),  and  also  In  (49),  f 
and  g  have  the  same  coefficient,  this  4X4  determinant  can  be  written  as  the 
sum  of  two  3X3  determinants  after  replacing  the  first  column  by  the  difference 
of  the  first  and  second  columns.  These  two  3X3  determinants  may  then  be  added 
to  form  a  single  3X3  determinant  after  absorbing  the  external  factors  Into  the 
top  rows.  Expanding  the  resulting  3X3  determinant  along  Its  top  row  and  setting 
the  result  equal  to  zero,  one  obtains  the  following  relation  between  the  fre- 
quency  and  the  wave  vector  k,  the  dispersion  relation  for  the  system: 
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In  order  to  put  ( 50 )  In  a  slightly  more  compact  form,  the  following  abbrevia¬ 
tions  are  Introduced. 


<•  5  dr 

poCv 


n  =  -  Index  of  refraction  of  unperturbed  medium 

oe  oe  r 

\  2 
n_  cE” 


I.  = 


oe  o 


L  -  — ^ — ■  =  power  In  Incident  laser  beam/unit  area 

e  -1  0 

A'  =  -2|—  AE2 
o  o 


B'  = 


c.  -1  E 


2^  i  (c2k2  .  eoeW2) 

5  (c2|t-  -  Eoe“-> 

The  dispersion  relation  can  now  be  expressed  In  the  form: 

j(u>  -  lK'k2)[w2  -  INwk2  -  (u2  -  A')k2]  -  (y-l)(u2-B' )u>k2 }  • 

r  -2  /,  .2  A  ,  2, 

,  e9/Mr 

Kj 


a  acne  a  1 
(a>  -  Itc'k2)  +  1  (u2  -  B')  + 

LAwk2+1^2.  iNwk2.  (u2.A,)k2jJ 

Eo  fc/.  +  tA) 

*  T  \  - ) 


Ak 
+  B 


poA 


(51) 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 


(58) 


For  comparison,  the  dispersion  relation  obtained  by  K.  A.  Brueckner  and 
S.  Jorna^  for  frequencies  outside  the  Raman  scattering  range  Is: 
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(59) 


u)[oj2  -  v2k2  -  INwk2]  CL  +  IL(vtu  +  16)  k2  =  0 


where 


v  <**  -  3“  W  • 


..*f£ 

n  J  c 
oe 

2Aev^ak? 
6  = — 


n‘c 

oe  p 


n  c 
oe 


4(  -?^)2  U 


—  k  ) 
n  *z' 
oe 


(59a) 


(59b) 


(59c) 


Although  equation  (58)  is  considerably  more  complicated  In  structure 
than  equation  (59),  the  general  features  of  the  two  equations  are  the  same. 
As  a  first  approach  to  the  analysis  of  (58),  one  should  realize  that  (53) 
relates  the  power  In  the  primary  laser  beam  to  E‘.  Thus,  the  free  modes  of 

*  '  i  ^ 

the  system  can  be  obtained  by  letting  E2  -*•  0  In  (58).  With  no  power  In  the 
Incident  beam,  therefore,  (58)  reduces  to 


(“'1^k2)w 
m  w 


2  -  INwk2  -  u2k2)  -  (y~l  )u2a>k2  * 


•  c2(»J  +  k2  +  a^kj)  -  el(0)  +  »>2* 

■ 

•  C^k2  +  k2  -  2kLk3)  -  eL(0)  (0)L  -  0))2  =  0  I 

The  first  factor  In  (60)  contains  a  non- propagating  thermal  wave  and  two 
damped  sound  waves  coupled  by  the  term  (y-l)uwk  .  The  last  two  factors 
correspond  to  the  four  free  modes  for  scattered,  undamped  electromagnetic 
waves: 


o+  %1  + 


where  a  =  +  1 
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2 

specified  by  given  values  of  k  and  k^.  Two  of  the  roots  are  low  frequency 
(u  «  jjl)  ,  whereas  the  other  two  have  frequencies  of  the  same  order  as  the 
laser  frequency.  It  Is  clear  that  the  roots  at  the  high  frequency  should  be 
eliminated,  because  it  has  been  assumed  previously.  In  evaluating  time  averages, 
that  the  perturbed  solutions  vary  much  more  slowly  than  the  optical  waves. 
Therefore,  factors  like 

cA2  -  («> 

will  be  replaced  by 

eL(0)[2uL3rcLk+  '  “+]  »  (63) 

where, 

cl  ~  TZZ  *  velocity  of  light  In  the  medium  (64) 

/ei.(o) 

Thus  the  free  modes  will  Include  three  thermal -sound  waves  and  two  electro¬ 
magnetic  waves. 

All  of  the  terms  In  ( 6(J  result  from  the  left  hand  side  of  (58)  because 
the  right  hand  side  Is  proportional  to  the  power  In  the  primary  laser  beam. 

Now,  as  the  power  In  the  lager  beam  Is  turned  on,  the  right  hand  side  couples 
the  five  free  modes  described  by  ( 60) .  Additional  tiny  coupling  arises  Inside 
the  left  hand  side  Itself  through  the  A1,  B',and  terms. 

For  detailed  consideration,  the  case  of  a  primary  laser  beam  at  10.6  u 
propagating  through  air  at  approximately  10°C  and  at  standard  pressure  will  be 
discussed.  The  numerical  values  for  the  parameters  appearing  In  the  dispersion 
relation  are:^73 

WL  =  1.773  x  1014  sec"1,  1^  =  5,920  cm"1  pQ  =  1.25  x  10"3  gm/cm3 
N  *  0.284  cm^/sec  8  =  3.67  x  10"3  deg"1 
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Cy  »  7.143  x  10  erg/gm-deg 


Y  *  1.4 


k'  ■  0.28  cnr/sec 


-4 


-4 


nfl  »  1  +  2.82  x  10 
u2  ■  8.39  x  108  (cm/sec)2 


atLe 

-  0.988  a 
p0Cv  o 


e  8  1  +  5.65  x  10 
o 

a  =  (3  x  ID"7  cm-1)  oo 

2 

an  cu  8  o 

—2 - -  (3.831  x  10*)  a 


2 

e2(0)Eo 


(1.37  x  10"25  -s~  j  )  I, 
gm  cm 


A  8  0.452  |£-  +  (1.10  x  10"22  )  IL 


gm 


%U  .  24.7  cgjdea  t  l6  w  „  10-21  cm  sec*  dea,  t 
5,1  gm 

eQe  8  1  +  5.65  x  10'4  +  (1.37  x  10"25  — )  IL 

gm  cm 

2  4 

A'  8  (2.84  x  10"15  -  )  IL  +  (1.38  x  10"36  ^“)  IL 

B  .  »  -  (4.84  x  10"28 - ^ - )  I, 

T0  gm  an  deg 


■o  9*11 

8*  8  -  (0.662  x  l(f36  ^  )  ij 


In  the  above  list  a  dimensionless  absorption  constant  aQ  of  order  unity  has 

been  Introduced  and  the  power  IL  Is  In  units  of  ergs/sec  per  square  centimeter. 

Now,  using  these  numerical  values,  one  finds  that  the  power  dependent  terms 

2 

are  very  small  for  power  fluxes  less  than  10  MW/cm  ,  except  for  the  term  which 
represents  energy  absorption.  In  other  places  In  the  dispersion  relation,  the 
power  dependent  terms  are  connected  with  the  nonlinear  Index,  and  will  be 
omitted  In  the  following.  (The  terms  omitted  are  related  to  self-focusing  In 
a  manner  described  by  Brueckner  and  Jornaf4^)  This  neglect  of  the  nonlinear 
Index  and  of  the  weak  dependence  of  the  optical  coefficients  of  gases  on  the 
temperature  for  fixed  density  allows  the  simplification  of  the  dispersion 
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relation  to: 


|[w  -  Itc'k2]^2  -  INwk2  -  u2k2]  -  (y-1)  u^k2}' 


•[cLk+-aj+][CLk_-co_]  « 


Ak  w,  I,  .  v 

’  “2^~  |[CLk+-w+]  +  CcLk_-a>_3|- 

2 

[prtA  ,  acn.Bu  n 

T  <*■*' *2> ♦ <  -r— J 


(65) 


poA 


In  addition,  for  air  at  reasonable  powers,  the  term  -S~  (w-lic'k  )  on  the  right 

acn^Bu2  c 


hand  side  Is  negligible  compared  to  — *= —  .  Introducing  the  variable 

Lv 

*4 


(66  ) 


Instead  of  k3,  and  defining 


T  = 


_ 


^cp 


(67  ) 


which  corresponds  to  the  parameter  used  In  Ref.  [4],  the  dispersion  relation 
can  be  written  In  the  form 


[(oMK'k2)(w2-1Nwk2-u2k2)  -  (y-1  )u2ujkZ]  3 


±LV2  k2 
2  vs  K 


(68) 


The  problem  at  this  stage  Is  the  determination  of  the  maximum  growth 
rate  of  any  Fourier  component  of  a  distortion  of  the  plane  wave  as  a  function 
of  the  absorbed  power  from  the  bean.  That  Is,  one  must  solve  the  dispersion 
relation  for  the  frequencies  as  a  function  of  k,  v,  IL,  and  the  characteristic 
parameters  of  the  medium,  and  then  find  the  maximum  value  of  (-  Im  o>)  for  real 
v  and  k  with  |v|  <1.  Such  a  problem  cannot  be  solved  analytically  without 
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further  approximations.  One  region  of  Interest  would  be  the  high  power  limit, 
where  the  driving  term  would  overwhelm  the  losses  resulting  from  thermal  con¬ 
duction,  viscosity,  and  the  absorption  of  electromagnetic  energy.  In  that 
case,  all  the  Imaginary  terms  In  equation  (68),  with  the  sole  exception  of 
the  "1"  Immediately  preceding  x,  can  be  dropped,  reducing  the  dispersion 
relation  to  the  form  employed  by  Brueckner  and  Jorna  In  equation  (45)  of  the 
reference  [3]. 

To  proceed  analytically,  Brueckner  and  Jorna^  neglected  the  second 
term  In  the  brackets  on  the  right  hand  side  of  (68),  and  assumed  that  the 
maximum  growth  rate  would  occur  somewhere  on  the  curve  In  the  v,  k  plane 
determined  by  the  constraint 


Along  that  curve,  the  maximum  growth  rate  Is 

("Im  w)max  (K08)  •  (70) 

which  corresponds  to 

vs k  -Vi[  osiM  VnW]  (71) 

(These  results  differ  from  those  In  reference  [4J,  which  are  erroneous.) 

There  Is  no  assurance  that  the  actual  maximum  growth  rate  does  lie  along  the 
one-dimensional  subset  of  the  v,k  plane  assumed  In  reference  [4].  In  reference 
[5),  a  search  was  conducted  along  the  line 
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Is  two  percent  smaller.  No  other  curve  In  the  v,  k  plane  has  been  found 
which  allows  an  analytical  search.  Nevertheless,  one  suspects  that  these 
answers  are  sufficiently  close  and  that  further  analytical  effort  Is  not 
justified,  because  of  the  previous  approximations. 

An  Interesting  unknown  not  discussed  previously  Is  the  power  flux 
required  to  stimulate  these  Instabilities.  This  threshold  power  Is  clearly  a 
critical  function  of  the  losses  In  the  system,  which  therefore  renders  It 
Important  to  treat  them  carefully.  If  the  second  term  on  the  right  hand  side 
of  equation  (68)  is  dropped,  the  Instability  appears  to  have  no  threshold, 
because  the  conduction  loss,  which  must  be  overcome,  vanishes  as  k  ■+•  0. 
However,  as  k  ■+  C,  the  Stokes  and  anti-Stokes  terms  on  the  right  hand  side  of 
equation  (68)  tend  to  cancel  each  other,  and,  therefore,  there  Is  a  threshold 
power  flux  for  these  stimulated  thermal  Rayleigh  scattering  Instabilities. 
This  threshold  was  determined  by  computer  to  be  312  ergs/sec  per  square 
centimeter.  (The  computer  program  Itself  Is  exhibited  In  Appendix  A  of  this 
report.) 

The  presence  of  a  wind  does  not  alter  the  growth  rates  for  distortions. 
This  can  be  easily  seen  by  considering  the  problem  from  a  frame  of  reference 
moving  with  the  fluid.  A  uniform  beam  remains  a  uniform  beam  In  the  moving 
frame,  although  Its  direction  of  propagation  Is  shifted.  This  shift  In 
direction  has  no  effect  upon  the  stability  discussion. 

In  order  to  attempt  to  examine  the  behavior  of  the  beam  In  more  detail, 
without  being  restricted  to  the  linearized  equations,  a  computer  solution  of 
the  full  set  of  laser-fluid  equations  was  attempted.  Certain  computational 
-.problems  were  encountered  and  are  discussed  In  the  next  section.  Then,  In 
Section  IV,  the  results  of  the  computer  calculation  are  discussed. 
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PART  III 


STABILITY  AND  UTILITY  ANALYSIS 


It  is  well  known  that  in  the  solution  of  partial 

differential  equations  two  problems  arise.  The  first 

problem  is  that  the  system  of  difference  equations  maj 

not  be  stable  against  error  growth.  The  second  problem 

is  that  the  solution  of  the  difference  equations  may  not 

converge  to  that  of  the  differential  equations.  Von 

Neumann  ^^roposed  a  necessary  condition  for  stability 

for  linear  partial  differential  equations  of  parabolic 

and  hyperbolic  types.  This  condition  has  been  gener¬ 
ic  [10] 

alized  by  Kata  and  others. 

In  this  section,  a  handy  sufficient  condition  for 
stability  of  difference  systems  associated  with  parabolic 
and  hyperbolic  equations,  both  linear  and  non-linear, 
will  be  discussed.  Consideration  will  be  given  to  the 
generation  of  errors  in  computer  solutions  and  a  new 

criterion  for  the  utility  of  numerical  integration 

•  . 

schemes  for  complicated  systems  of  equations  will  be 
proposed.  The  new  term  ‘’utility"  is  introduced  to 
distinguish  the  criterion  from  the  better  known 
criterion  of  "stability".  Some  simple  applications  of 
this  criterion  will  be  examined  in  this  chapter,  as  a 


prelude  to  its  use  in  the  more  complicated  laser-fluid 
problem. 

The  concept  of  “utility"  is  basically  very  simple 
A  differencing  scheme  is  "useful"  if  the  solutions  of 
the  difference  equations  are  sufficiently  accurate 
approximations  to  the  solutions  of  the  partial  difl’eren 
tial  equations  over  the  domain  of  interest.  Such  an 
integration  technique  may  not  be  stable,  according  to 
the  usual  definition  of  "stability"  and,  thus,  may  not 
be  convergent  either.  This  lack  of  stability  and 
convergence  does  not  cause  alarm,  because  the  accuracy 
can  be  estimated.  The  best  analogy  would  be  an 
asymptotic  series  expansion  for  a  function.  Such  an 
expansion  may  not  converge  to  the  function,  yet  can 
still  be  a  highly  satisfactory  description  of  the 
function  over  the  domain  of  interest. 


I 

A,  The  Stability  Criteria 

The  most  general  form  for  a  multi-level  difference 
scheme  is 

+  + . +  A,W„+fc „=0  (73) 

where  - A,  are  finite-difference  operators  which 

are  functions  of  time  in  general.  The  time  is  designated 
by  Yl  and  Uf*  is  the  dependent  variable  of  the  difference 
equation.  The  actual  solutions,  W*  sssu^ma*)  $  for  the 
associated  system  of  differential  equations  will  not  be 
the  same  as  because  of  truncation  errors  h*  introduced 
at  each  step  of  the  calculation. 

An  equation  of  the  form  shown  in  ( 73)  results  from 
converting  any  system  of  partial  differential  equations 
into  a  system  of  difference  equations.  The  vector 
is  a  vector  with  many  components  representing  the  values 
of  each  of  the  dependent  variables  at  each  of  the 
spatial  mesh  points.  In  addition,  some  of  the  components 
of  can  be  taken  to  represent  temporal  derivatives 
of  the  dependent  variables  if  the  equation  is  of 
second  or  higher  order  in  time. 

Assuming  that  (  73) ,  when  supplemented  by  suitable 
boundary  and  initial  condition  has  a  unique  solution, 
one  may  write 
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—  ^-iu£tj-r  +  - . +  +  J«  ,  (74) 

whore 

Bo  ■= 

The  solution  of  the  system  of  differential  equa¬ 
tions,  evaluated  at  the  time  r\4t  ,  will  be  denoted  by 

,  ard  the  value  of  the  solutions  of  the  difference 
equations  generated  by  the  computer  by  Uf*  •  The  differ¬ 
ence  between  these  values  will  be  called  the  error  and 
will  be  denoted  by  •  In  general  because  of 

round  off  errors. 

If  £„  is  sufficiently  small  then  it  will  satisfy 
a  linear  system  of  difference  equations. 

e**l  =  C%.\  + . +  C,eh  t  <j„  .  (  75) 

(xn  represents  the  local  errors  introduced  at  each  step 
and  is  the  sum  of  truncation  errors  and  round  off  errors. 
For  linear  equations,  and  will  be  independent 

of  M  if  the  time  does  not  enter  explicitly  into  the 
linear  equation.  However,  if  the  basic  system  of 
equations  is  non-linear,  then  C;(k)  &C*0  . 

The  vectors  . ,  as  well  as 

W*» ,  Uf„tl . tuXtg-ican  be  arranged  so  that  the  com- 
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ponents  form  new  vectors,  £*  and  W„  ,  with  g-  times  as 
many  components 


/ 

I  • 

and  W*  — *• 

Wi«-  i-i 
• 

• 

\Ci  ) 

_ 

1  Vn  y 

(76) 


Then  equations  (  74)  and  (75)  can  be  rewritten  in  the 
form 


BnW«  +  J 


X  , 


(77) 


En-*i 

where 


C*  “h  (sfr 


sss 


cl  ss  CC"/0*)  = 


B^*) - 8,(i)  £  J 
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o 
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o 

t 

« 

I 
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o 

I 


ch(^ 

0 

I 

I  V 

% 

I 

I 
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o 
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o 
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c,(*)  cy) 

o  o 


-  0 


-  1 


0 

1 

1 

l 

1 

I 


0  J 


(78) 


(79) 


(80) 


and 


The  matrices,  I,  appearing  in  (79)  and  (80)  are 
identity  operators  for  the  subspaces  over  which  the 
difference  operators  13^  and  C:  operate. 

For  homogeneous  difference  equations  with  time 
independent  coefficients,  the  index  71  on  o„can  be  dropped 
and  equation  (77)  becomes 

W„t,  **  6Wn  (82) 

where  6  •  Kant  orowitch^11  defined  a  system  of 

linear  homogeneous  difference  equations  with  time 
independent  coefficients  to  be  "stable"  if,  for  specified 

/v  ■  M 

T,  there  exists  a  positive  number 't  such  that  L 
is  uniformly  bounded  for  all  such  that  omk  t'Y 
and  for  all  integers  n  such  that  0***fc<  T  , 

A  straightforward  generalization  for  the  case 

•  A/ 

where  B  is  also  a  function  of  the  timo  index  n  would 
bo  the  following. 

Definition;  A  system  of  difference  equations  is  colled 
stable  if,  for  specified  T,  there  exists  a  positive 
number  T  such  that  the  set  of  operator  products 

S(h'lA*) - . is  uniformly  bounded  for  all 
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£y ^  such  that  o<6*4.'t  and  for  all  integers  vt  such 
that  o<m<i i<T  . 

In  case  the  spatial  coordinates  do  not  enter 
explicitly,  tht  analyses  of  the  system  can  be  simplified 
by  a  Fourier  transform  technique.  In  such  cases  equation 
(82)  will  lead  to  an  equation  of  the  form 

—  $(**,£)  (83) 

where  is  the  spatial  Fourier  transform  of  W*  . 

The  matrix  dr is  called  the  amplification  matrix. 

An  equivalent  statement  of  the  Kantorowitch  stability 
criterion  is  that  for  specified  T,  there  exists  a  posi¬ 
tive  number  c  such  that  the  set  of  operators  q'cat,#) 
is  uniformly  bounded  for  all  values  of  and  for  all 
such  that  and  for  all  integers  n  such  that 

o  c.  haJf  c  T 

The  value  of  the  Kantorowitch  stability  criterion 
lies  in  the  proof  by  P.D.  LcuP0,^that  in  the  limit 
that  tends  to  zero,  all  systems  which  are  stable  lead 
to  solutions  which  converge  to  those  of  the  associated 
differential  equations  and  vico  versa.  To  proceed  it  is 
useful  to  introduce  some  standard  terms: 

(1)  The  Lp  norm  of  a  column  vector  \f  ,  denoted  by  jvL 

T  J6  * 

is  defined  by  I  Vlj>  ^  jtfcfj  • 

»  • 

(2)  The  maximum  norm  or  the  norm  of  a  column  vector. 


V  •,  denoted  by  is  defined  by  I  U'tp  **  I  ^  I  • 

(3)  The  induced  matrix  norm  of  a  matrix  ,  denoted  by 

Iffcilp  or  l|  is  defined  by  I  £||  *■  ^77^ 

(4)  The  radius  of  a  matrix  £  ,  denoted  by  ft(.X) ,  is 
defined  by  RtjOss  ***(  %  |  *^| .)  — 

The  radius  of  a  product  of  matrices  is  not  greater 
than  the  product  of  the  radii  of  the  matrices, 

(5)  The  spectral  radius  of  a  matrix  ^  ,  denoted  by 

f(*)  ,  is  defined  by 

whore  (a**  1,2,  •■•,**)  are  the  eigenvalues  of 
J6  •  Clearly,  it  is  true  that  f(£)^RC£) 

If  X  is  a  normal  matrix  (commutes  with  its 
hermitoan  adjoint),  then  it  oan  easily  be  shown 
that 

(fftjf-ftcrj-ut*)]’! 

<• 

With  the  aid  of  these  terms,  the  Von  Neumann 
necessary  condition  for  stability  can  be  stated  in  the 
following  way.  A  necessary  condition  for  stability  is 
that  there  exists  positive  numbers  'V  and  D  such  that 
the  spectral  radius  of  the  amplification  matrix 

■&)  satisfies  I  for  all  and  all 

-oir  such  that  ,  if  is  a  noi'mal 

matrix  (commutes  with  its  hermition  conjugate),  then  the 
Von  Neumann  condition  is  sufficient  as  well  as  necessary. 
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Many  sufficient  conditions  for  stability  have  been 
proposed 12^  The  analysis  for  stability  has  been 
generalized  to  the  case  in  which  the  coefficients  of  the 
differential  equations  are  time  dependent  by  P*D.  LJ10.13] 
and  many  others  E10-!  It  was  further  generalized  by  W.G. 
Strang ^  14lnd  many  othercf^io  the  case  in  which  the 
fundamental  system  of  equations  is  non-linear.  These 
extensions  will  not  be  detailed  here.  Instead  a  handy 
sufficient  condition  for  stability  will  be  considered. 
Theorem:  A  sufficient  condition  for  stability  is  that, 
for  specified  T,  there  exist  positive  numbers  rand  D 
such  that 

fc  ±  l-t  (84  ) 

for  all  Lt  such  that  ,  whore 

r  -  ***  a  ( S  ■*  I  (85) 

for  all  positive  Integers  n  <■  T/&t.  For  a  linear 
equation  with  time  independent  coefficients, 

R=RCS)«ligC. 

Proof:  For  R  *  J  ■+■  U>*t  ,  it  follows  that  R " ^ 
where  K  » epr  .  Therefore, 

II  *5(11-1  ,  aX-)  •  *  •  * 

So,**)) 
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^  R(8(n,4*))  - R(&  a  A*)) 

for  all  At  and  M  such  that  OA&tt't  and  OAnaJtA  T 
so  that  the  maximum  norm  is  uniformly  bounded*  Since 
the  maximum  norm  is  uniformly  bounded,  then  the  matrix 
product  is  uniformly  bounded.  Thus  the  stability 
definition  given  above  is  satisfied,  Q.E.D. 

A  matrix  is  called  a  dominant  diagonal  matrix 
if  the  "interior  radius",  R-,n  (V)  >  0  ,  where 

■  mT  I  KJ  "  M«*  i  (86) 

Lemma:  For  a  two  level  implicit  scheme  /U**tO i/,,+, «=  f 
a  sufficient  condition  for  stability  for  specified  T 
is  that  there  exist  positive  numbers  X  and  Q  such  that 

0  <  -Jr  ^  I  +  ,  C  87  ) 

Ki*  * 

for  all  such  that  0  <  T  and  for  all  integers 

n  such  that  o  <  <  T  »  v/here 


Proof;  Since  f?ir» (AW)>0  for  all  y\  such  that 

it  is  true  that  B  (  I  +  G  +  El  +  £* ■*•••*•  )  whore 


X  o 


...  \ 


8 


o  A! 


e  =  i-/»b  ; 


;  t«;/  eif  /  0£j 
(. «« 


•for 

for  ?»,). 


It  follows  that 


«(•*%>)£  R(6)[  I  +  R(£)  +  (R(B»'  +  CRCS)j+  -■■■]—  $Sij 

3  .  ■  2 


i«sj 


_ I. 


Therefore,  since  —  K  inhere  K  =  <^r 

it  foil  03/s  that 

II A ",0)A''(t)--A%-'M‘h)ile0 

*  '  . 

•  —  Vo/f'oo) 

t&  R(A"(»)R(A i'M)  -  RCA"(»-'»Rl/f(n)) 

4  _J _ I_ . . 2 _ L_ 

RUAO)  R;JA(i))  1Uh("-o)  RfA(»)) 


I 


Thus  a  uniform  bound  is  obtained  and  the  stability 
definition  is  satisfied.  Q.E.D. 

The  above  lemma  offers  a  sufficient  condition  for 
stability  of  an  implicit  difference  scheme.  It  is 
conceivable,  of  course,  that  some  stable  implicit  differ¬ 
ence  schemes  might  not  satisfy  the  criterion  given  in 
this  lemma.  More  generous  sufficient  conditions  can  be 
given,  but  none  is  as  easy  to  apply  as  this  one. 

As  a  practical  note,  it  might  be  remarked  that, 
for  actual  computation,  one  is  often  interested  in  using 
the  largest  value  of  Ak  consistent  with  the  stability 
requirements  for  a  given  mesh  size.  This  is  certainly 
the  situation  in  the  computer  solution  of  the  laser»fluid 
equations  attempted  in  Part  IV  of  this  study.  In 
such  cases  the  relation  of  X  to  the  spatial  mesh  size 
is  of  great  interest.  The  spatial  intervals  are  often 
dictated  by  the  scale  of  the  initial  configuration. 

For  example,  if  the  dependent  variable  at  t-0  is&~'z,fZ 
then  it  is  unlikely  that  values  of  AX  larger  then  0.3 
would  be  useful  from  the  standpoint  of  convergence 
requirements.  The  connection  between  X  and  the  spatial 
mesh  size  can  be  illustrated  with  the  diffusion  equation 


where  the  term  -4- it  is  included  to  allow  some  exponential 
growth  in  the  sclution.  The  simplest  explicit  difference 
*  scheme  representing  (  88)  is 


-2U*  + 

d 


ui-<] > 


(89  ) 


■  V  . 

where  n  is  the  time  index  and  is  the  spatial  step 
index.  This  difference  equation  can  be  put  in  the  form 


(W 


+ui] 


(  90) 


Using,  for  simplicity,  the  handy  sufficient 
criterion  given  by  the  theorem  above,  one  would  use  ( 90) 
to  evaluate  k  •  Since  the  time  does  not  enter  explicit¬ 
ly,  R  =  R.  •  It  is  trivial  to  compute  R.  directly  from 
(90)  without  writing  down  the  matrix  8  explicitly. 

Doing  so,  one  finds 


*-|i- ^  +  **1  +  1^1  + 1^1. 

Since  K  ,  Ajt  ,  andC^X)  are  all  positive,  this  expression 
simplifies  immediately  to 


The  stability  condition  now  enters  through  the  argument 
necessary  to  remove  the  magnitude  bars  in  (91  ) .  One 
must  be  careful  to  realise  that  the  D  in  the  stability 
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theorem:  £  /  +  P4f must  not  be  a  function  of  ,  All 

of  stability  analysis  is  couched  in  a  framework  in  which 

AX  is  supposed  to  be  allowed  to  go  to  zero  in  some 

manner.  The  D  of  the  theorem  must  not  be  affected  by 

Ax  —>  0  .  In  other  words,  one  is  not  allowed  to  consider 

■  ^r~rnr  to  be  of  order  Air  as  it  would  be  for  fixed  . 

t-ax) 

The  theorem  does  not  intend  ax  to  be  fixed.  Thus  the 
removal  of  the  magnitude  bars  in  (91  )  requires  further 
analysis.  The  term  6  may,  of  course,  be  considered 

to  be  of  order  a*  because  'O-  is  indeed  fixed,  independent 
of  ax  .  If  it  is  true  that 


2  7<a£  y  i 

(axX  ~  •* 


then  (  91 )  becomes 


K  t  1  -  —  /  -f 


(92) 


(93) 


and  plays  the  role  of  V  in  the  theorem  and  the  scheme 
is  stable.  The  stability  condition  is  (  92)  ,  so  that 

At  £  53 't  (94) 


Thus  7T  has  been  discovered  and  related  to  AX  ,  There 
is  one  other  possibility  for  removing  the  magnitude 
bars  in  (91  ) .  If  -——jr  >  2  ,  then  ( 91  )  becomes 


so  that  the  scheme  fails  to  meet  the  sufficient  con- 

« 

dition  given  in  the  theorem. 

The  points  to  be  emphasized  in  the  example  above 
are  the  following.  The  quantity  D  appearing  in  the 
stability  theorems  must  not  depend  on  the  spatial  mesh 
sizes.  The  usefulness  of  the  sufficient  condition 
arises  through  the  relation  derived  connecting  't  with 
the  spatial  mesh  sizes.  This  relationship  typically 
arises  by  requiring  the  part  of  R  that  is  not  of  order 
Air  to  be  less  than  or  equal  to  one.  Furthermore,  a 
term  containing  spatial  mesh  sizes  in  the  denominator 
is  not  of  order  &1:  .  As  indicated  earlier,  there  are 
certainly  more  generous  sufficient  conditions  available 
in  the  literature  than  the  one  illustrated  here,  but 
none  are  as  quick  and  easy  to  apply.  Since,  however, 
in  computational  work  one  never  actually  lets  AX  and 
go  to  zero,  it  is  not  clear  that  analysis  based  on  such 
limit  notions  is  necessarily  the  best  approach  to 
practical  problems.  As  a  matter  of  fact,  peculiar 
paradoxes  such  as  the  one  discussed  on  page  230  of 
Ref erence[10], arise  from  such  limit  notions.  In  hope  of 
obtaining  criteria  of  greater  practical  value,  the 
utility  viewpoint  has  been  devised  and  will  be  presented 
next. 
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B.  The  Utility  Criteria 


The  theorems  mentioned  above  are  not  practically 
applicable  in  general  because,  by  omitting  local  error, 
there  is  no  way  to  estimate  the  maximum  error  growth. 
Also,  the  increment  of  space  variables 
and  the  increment  of  time  variable  do  not  tend  to 
zero  in  actual  computation.  Therefore,  a  new  criterion, 
"utility"  is  proposed  here  to  avoid  some  unnecessarily 
stringent  requirements  and  to  limit  the  maximum  error 
growth  at  the  same  time.  The  notation  of  the  preceeding 
section  will  be  employed  below. 

timv  For  given  6  and  T  ,  and  for  any  given 
initial  and  boundary  conditions,  a  difference  scheme  is 
useful  if  there  exist  positive  numbers  V,  and  %  such 
that  for  any  a k  in  the  interval  X&ajr&Z 'L  t  the 
absolute  errors  are  bounded  by  6  ,  for  all  integers  m 
such  that,  o  f  where  T  is  the  physical  time  of 

interest  for  a  given  problem,  is  the  increment  in 
time,  ki  is  the  time  index,  and  T;  and  Zx  are  to  be 
determined  by  the  round  off  and  truncation  error, 
respectively. 

Theorem;  For  given  6  and  T  ,  a  necessary  and  suffi¬ 
cient  condition  for  utility  is  that  there  exist  positive 
numbers  T;  and  T*.  such  that  for  any  -a  k  in  the  interval 
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'Ci  6  i  'ti 

c*{  C*,.,  [  C, ,.,(  ••■•  (  C,£,  +4»)  +  "’•  H£„.J +  *„.,}  +  *,][ 

-  *  » 


(  95) 


for  all  integers  n  such  that  o  j  -f  where 

\ZlS„l----]+Z~}+lL 

=  ^ |  [  •**•]  t +£wj  #  ( 96) 

Proof :  From  the  difference  equation  for  the  error,  one 
gets 

El~C,El  +  Gll 

t  C  i  {  C,  b,  t  }  +  *^ 


E4.  +  *; 


Eh*,5—  ?*{ 5,.,  [••••(  CE,  +  *,  )---*  }  +  *n  • 


It  is  easily  seen  from  the  definition  of  £„  that 


Max 

1* 


Kl  =  Mr  i£»L, , 


(  97) 


for  all  integers  n  such  that  o4m44T  •  Therefore, 
6  is  the  upper  hound  of  the  maximum  absolute  error 
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The  maxima 


for  any  in  the  interval  V  . 

shown  in  (97  )  exist  because  only  a  finite  set  of  values 
.of  n  are  involved  for  t;  ?0  ,  Q.E.D. 

For  linear  equations  with  constant  coefficients, 

Z  is  independent  of  n.  If  one  makes  the  additional 
assumption  that  <5r  is  also  independent  of  n,  then 
equation  (95)  becomes 

lei  +c  +  c*+  ^5+  evzl, 


/vMI 

c  - 

I 

c  — 

X 

The  local  error  <3rn  can  be  expressed  in  the  form: 

S-H 


(?rn  —  pn  +  H„  *  4  At. 


(98) 


where  ft  and  H«  are  the  local  round  off  error  and 
truncation  error,  respectively.  For  a  general  differ¬ 
ence  scheme,  one  may  write  H*  =h«aiul ,  where  s  is 
determined  by  the  difference  scheme.^ 

A  handy  sufficient  condition  for  utility  will  now 
be  given.  The  following  notation  is  introduced  for 
convenience: 


where  the 'maxima  are  computed  over  the  finite  set  of 
values  of  n  such  that  o  <.  for  specified  7  ,  T,  , 

and  Tt  such  that  %  -  'ti.  *  The  set  of  integers  n 

is  finite  because  of  the  upper  bound:  -J~  • 

Theorem:  For  given  6  and  J  ,  a  sufficient  condition 
for  utility  is  that  there  exist  positive  numbers  T, 
and  Tt  such  that 
S’"*1  / 

-  (  p  (102) 

for  all  Ajt  such  that  %  £^*£ti'and  for  all  integers 
n  such  that  o*  T  , 

Proof :  The  quantity  appearing  in  line  (  95)  of  the 
preceeding  theorem  is 

I~  r  ft.  -  ft,  A/  ft#  ft  #V  I  #V  | 

c"{C,-,f  C*  (■—  +4,  )  +  ----)+^nwJLJ  +^n-i  f 

✓ 

*a=:R.{R(R( -  (R^+(t)h - 

N  jtlMtV 

=  £<  +R  <«‘4?+ - +R"j? 

—"♦I  . 

^-j~-(P+un 

Q.E. D. 
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For-  vt  =  ,  the  LHS  of  (102)  becomes 

j  - — C  p  +  Ji&t**')  .  If  is  very  large,  then  the 
truncation  error  IT  will  become  large  and  the 

error  can  surpass  the  bound.  If,  on  the  other  hand,  Afc 
is  very  small,  then  the  round  off  error  from  ~P  can 
build  to  a  large  value  and  surpass  the  error  bound  6 r  . 
There  may,  however,  exist  a  region  of  a**,  C  % tl )  , 

in  which  the  utility  condition  is  satisfied.  Obviously, 
if  Rs  |  ,  then  the  proof  above  breaks  down  in  the  last 
line  and  the  criterion  can  be  stated 

(yi+O  c  P  +  ♦  (103) 

for  all  n  such  that  .  This  criterion  would 

require  that  the  minimum  of  -Jy-  +  be  less  than  k/ r  . 
The  condition,  (103  )  would  not  bo  useful  unless  one 
actually  could  estimate  p  and  k  .  This  can  be  done 
easily  only  for  very  simple  schemes.  There  is,  however, 
no  necessity  to  take  R-|  since,  without  loss  of  gener¬ 
ality,  a  larger  value  can  be  used  in  place  of  R  •  Thus 
(102)  may  be  used  in  all  cases, 

A  very  useful  corollary  to  this  theorem  can  be 
given,  based  on  the  precondition  that  ||  (T  Ij^  can  be 
written  in  the  form 

||  1  +  tTh^  so  that  R=  I  +  (104  ) 
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where,  unlike  the  case  of  stability  criteria, 

(  and  hence  $  ).  is  allowed  to  depend  on  the  spatial 

mesh  sizes.  One  should  notice  that  condition  (104  )  for 
R  is  not  the  same  as  the  stability  criterion  (86  )  for 
two  reasons:  $  ,  unlike  D,  is  allowed  to  depend  on  the 
spatial  mesh  sizes  and,  furthermore,  the  R  in  (104  )  is 
equal  to  II  C*  only  for  ,  a  finite  set  of 

values. 

Corollary:  If  precondition  (  104)  is  satisfied,  then  the 
sufficient  condition  (102)  becomes 

(LLSzV-J  (  P  t  W*1;  —  6* .  (105) 

$  &'k 


The  precondition  bears  further  discussion.  Certain¬ 
ly,  for  the  difference  scheme 


26?  -h  £4*. 
i  GK 


(106) 


it  is  true  that  l|  C'|[o==2 ,  so  that  the  precondition 
( 100)  is  not  satisfied.  In  such  a  situation,  the  suf¬ 
ficient  condition  (105  )  could  not  be  applied  and  one 
would  have  to  use  (102),  Equation  ( 106)  represents  an 
unusual  situation,  however.  For  difference  equations 
resulting  from  differential  equations  written  in  a 
form  where  only  first  derivatives  with  respect  to  time 
arise,  then  one  can  have  the  explicit  scheme: 
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(107) 


£**'  -  E" 


where  Qn is  a  matrix  operating  on  the  components  of  fcw  , 
Notice  that  in  this  cftnonical  situation  an  unusual  factor 
like  the  2  in  (106  )  does  not  enter. 

It  happens  that  in  very  many  difference  schemes, 
the  precondition  ( 1 04 )  is  automatically  satisfied. 

There  are  cases,  however,  for  which  (104  )  imposes  a 
real  condition.  For  example,  the  diffusion  equation 
(88)  leads  to  (91): 


Thus  I i  £*(1^  can  be  put  in  the  form  /  +  as  shown 

in  (93)  with  Z  =-£,  only  at  the  price  of  accepting 
the  condition 


2  fra*  *  i  (108) 

( ax )l  " 

which  happens  to  be  the  stability  condition  (92).  Thus, 
in  this  example,  one  gets  the  stability  condition  as  a 
precondition  for  the  utility  condition; 

— ( p+v^1)  (lo9) 

Ar &"K 

and  furthermore  $  is  independent  of  mesh  aiae  in  this 
case.  Roth  (108  )  and  (109  )  must  be  satisfied  in  order 
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to  appeal  .to  the  sufficiency  condition  (  105)  for  utility. 
Instead,  one  can,  of  course,  use  the  simple  condition 
(102)  If  •£*=0  in  (88)  ,  then  one  still  gets  the 
stability  criterion  as  a  precondition  and  one  can  use 
(  103)  for  the  utility  criterion  as  already  discussed. 

There  is  a  third  option  available,  if  one  wishes 
to  use  neither  (l 02)  nor  the  combination  ( 1 04  )  —  C 1 05  ) . 

One  can  decide  to  accept  a  slightly  more  restrictive, 
condition.  For  example,  from  (91)  one  can  write 


K  —  1  + 1  u>x)'  +  4] —  1  +% ft &i:  (no  ) 

and  use  Feff  in  (105)  to  get  the  rather  strong  condition 
for  utility 


im 

~  l,  (  p  + 


<  111  ) 


The  Schrodinger  equation  is  similar  to  the 

diffusion  equation  (88),  except  that  the  factor,  i, 

multiplies  ^  on  the  left  hand  side  of  the  equation. 

This  innocuous  looking  factor  of  i  has  the  crucial 

effect  that  no  explicit  difference  scheme  for  the 

[151 

Schrodinger  equation  can  be  stable.  This  can  easily 
be  seen  by  computing  R  .  This  time  there  is  no 
stability  condition  such  as  (92)  and,  correspondingly, 
the  precondition  C104I  is,  more  or  less,  automatically 


if 


satisfied  and  (  105)  gives  a  utility  criterion  which, 
unlike  (109),  involves  the  spatial  mesh  sizes.  The 
Schrodinger  equation  and  other  examples  will  be  analyzed 
in  detail  at  the  end  of  this  section  to  further  illus¬ 
trate  the  application  of  the  utility  criterion. 

The  practical  value  of  criterion  (105)  or  (111  ) 

lies  in  the  fact  that  large  values  of  n  will  often  be 
used,  say  n  £  A/  It  is  then  useful  to  rev/rite 

(105)  in  the  form 


(  I  +  I  +  (  e 


C  S**)  — 


(p  -f  W”)  <R4Wf')  ' 

where  it  is  presumed  that  the  error  bound  is 

generous  enough  compared  to  the  error."  in  any  single 


(112) 


step  that  the  1  will  be  insignificant  and  can  be  dropped. 
Taking  logarithms,  one  then  obtains  the  critical  condi¬ 
tion 


Cm  (I  +  hi  <*$£>.  +  — i—  Ujf\ - S_1 

a  /V  +  I  N+\  H  p  +  TyAl.**1  J 


~  [ 


'J. 


(113) 


p  -t  Ww 

where,  again,  it  is  presumed  that  the  error  bound  is 
generous  enough  that  a  term  can  be  dropped  on  the  8HS, 
It  is  also  presumed  that  hJ»\  so  that  N+l  •  Thus 

the  useful  form  of  (  105)  under  these  conditions  is 

(1,4> 
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Now  attention  is  called  to  the  fact  that  it  does  not 


matter  very  much  what  one  uses  in  (114)  for  the  ratio, 

. p  +  ,  because  only  times  the  logarithm 

of  this  ratio  enters.  For  example,  if 


loo 


and 


A/—  100, 


p  •+  hAi**' 

where  A/  is  the  number  of  steps  to  be  used,  then 
Aoq  Cl  +  «.  O.  a4  l>o(, . 

Then,  using  a  table  of  natural  logarithms,  one  finds: 


2.1  ar 


(115) 


If,  on  the  other  hand,  a  considerably  more  generous 
bound  is  used,  say, 


it 


P  + 


At  10 


for 


bj  = 


•then 


A-JO  +  *=  O.tro  l 


so  that  the  condition  is 


I 


**  *  (f)T  • 


20 


Thus  changing  the  error  bound,  b  ,  by  a  factor  of  10 
only  relaxes  the  criterion  for  by  a  factor  of  14. 
Thus  one  might  expect  to  try  to  use  the  rough  criterion 


-rr  with  Zo^to  , 


(116) 


for  almost  any  system  of  equations.  If  (116)  leads  to 
instabilities  in  the  computer  solution,  then  one  might 
increase  to  to,  say,  ho  and  try  again.  If,  on  the 
other  hand  (  116)  leads  to  stable  solutions,  but  is  too 
re.-trictive,  then  one  might  decrease  to,  say,  1  and  try 
again  to  get  stable  solutions.  Following  this  course, 

P  and  k  would  have  to  be  estimated  only  if  one  wished 
to  know  the  value  or  the  error  bound  £•  corresponding  to 
(1>6).  One  must,  of  course,  be  able  to  distinguish 

bo tweeil  iiic»l/iioiiiaLiet»l  iuStabll i  and  true  ]>hy  W  .i  Otll 

instabilities  (such  as  the  e  in  (88)  )which  may  arise 

in  the  computer  output.  With  some  experience  one  is 

able  to  spot  earmarks  of  some  of  the  mathematical 

[10] 

instabilities  quite  easily. 

For  implicit  two  level  schemes,  it  is  easier  to 
use  the  slightly  different  criterion  given  in  the 
following  lemma. 

Lemma:  For  a  two  level  implicit  scheme  A <"t!)EUtc=r  E*  ,  a 
sufficient  condition  for  utility  is 


0  < 


lu - T-fp  +  A^O^c- , 


(117) 
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for  all  integers  n  such  that  T  ,  where 

R;„<A)  es  wij»|  Aaa|  -  (117a) 

and 

R;,,*3  ^  (117b) 

In  some  ways  the  utility  criterion  is  similar  to 
the  "practical  stability  condition"  discussed  in 
Reference  but  the  philosophy  is  drastically  different 
due  to  the  fact  that  the  utility  criterion  makes  no 
attempt  to  deal  overtly  with  questions  of  stability. 

It  should  also  be  mentioned  that  the  utility  criterion 
is  easy  to  apply,  even  to  the  complicated  set  of 
equations  describing  the  laser-fluid  system.  The 
"practical  stability  criterion"  and  other  common 
stability  criteria  are  extremely  difficult  to  apply. 

Partly  the  difference  in  the  ease  of  application 
of. the  various  tests  arises  from  the  fact  that  stability 
and  bounded  error  growth  (utility)  really  are  slightly 
different  concepts.  In  large  part,  however,  the  dif¬ 
ferences  arise  because  some  of  the  stability  methods 
are  made  very  complicated  in  order  to  weaken  the  con¬ 
dition  on  &£  as  much  as  possible.  In  this  regard  it  is 
quickly  admitted  that  more  complicated  utility  tests 
*  can  easily  be  devised,  but  such  matters  will  not  be 


-55- 


considered  here.  It  is  a  part  of  the  utility  philosophy 
that  utility  criteria  should  be  less  complicated  than 
the  equations  one  wishes  to  solve  in  the  first  place. 

In  this  way  it  is  hoped  that  the  major  part  of  the 
computer  time  can  be  spent  solving  the  equations  of 
interest,  rather  than  trying  to  unravel  stability  or 
utility  criteria.  Stable  or  unstable,  convergent  or 
divergent,  if  the  computer  solution  correctly  describes 
the  phenomena  for  which  the  differential,  equations  were 
written,  then  that  is  utility  at  its  finest. 
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C.  Examples  of  Utility  Analysis 

In  order  to  Illustrate  the  application  of  the 
utility  criterion,  several  examples  will  now  be  given. 

The  criterion  will  be  applied  to  the  Schrodinger  equa¬ 
tion,  the  Dirac  equation,  and  a  nonlinear  partial  dif¬ 
ferential  equation  similar  to  the  laser  self-focusing 
equation.  In  each  case,  three  different  iteration 
schemes  are  considered:  a  two  level  explicit  differ¬ 
ence  scheme,  a  two  level  implicit  difference  scheme,  and 
a  simple  multi-level  '’Leap-Frog11  scheme.  The  method  of 
finite  differences  is  not  necessarily  the  best  means  of 
solving  these  particular  differential  equations,  nor 
are  the  selected  difference  schemes  necessarily  the 
best.  These  examples  are  presented  chiefly  to  illustrate 
the  application  of  utility  criteria.  For  simplicity,  it 
is  assumed  in  these  examples  that  £  -  /O  h  and  P  -  10 
Certain  transformations  are  introduced  in  the 
examples  in  order  to  reduce  truncation  errors.  This 
technique  is  a  standard  numerical  procedure  and  is  often 
called  the  "integral  method."  The  Schrodinger  equation, 
for  example,  may  be  put  in  the  form: 

&[yef  (-ife  Cv<Lt')j 

=  at  (~  A  fat ')  (ns) 
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Integrating,  now,  from  to  ,  one  obtains: 

--ikC'vdt) 

,  -ton 

=  I  vM)  +4  J  *H?W*p(-k 

Taking  if  e*  tf  &%p  (  JL  and  making 

1  An  “{j 

the  approximation  that 
tffVD  ^  x  p  c  - 

«  (VV")  c*P(-4  fVrfi'J, 

one  can  proceed  to  obtain  difference  equations  for  the 
better  dependent  variable,  .  The  integral  method  is 

incorporated  into  the  examples  be  lev:  for  the  Schrodinger 
equation. 

1.  Schroedlnger  Equation 


[vtf  J  019) 


i.e. 

H  ^  +  V/-  . 

Let 

if  sob  (^,  <2Xp  (  ^  J  V/r-M'  ) 


(120) 


(121  ) 


(122) 
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so  that 


t 


(  123) 


(a).  Two  level  explicit  difference  scheme  (forward  time 
differences ) 

4  %i*  **P<7kft  *  -  vs-,*)**)  -2 

+  f  e*P  ( A  [ ivtjtH  - 1 (ijt )  it) 

+  ?;  “H  e*P(-^f<  ^ -2 024 ) 

* 

Since  this  is  a  linear  equation,  the  error  satisfies 
the  relation: 

£ r*  ^  •  L  p  ^ 

4E*<«exP0*  ((hj4~  'tjj.')'**)- 

+1  t  **?<■£  ft  fy4-\ti4)  it) 
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( 125) 


+  f  A*  6  *^75  [  <  Vu  «*'  ~  ) 

+  p".  +  tt" 

*  d  «  ■*  j  K  , 


For  given  6  and  T  ,  the  necessary  and  sufficient  con¬ 
dition  for  utility,  involving  matrix  multiplication  of 
the  Cm  ,  can  be  determined  by  using  a  computer.  However, 
the  sufficient  condition  005  )  for  utility  of  AJtr  for 
fixed  N  =>foc>  where  T -  N&Jt**toc*Xcan  be  obtained  very 
easily  as  follows: 


RCn-~  I  +  Ut  —  |  I  +  + 

[  i^K)1  +  +  a/ 1 


^ +  ^  (  (ax/  +  T^y?  + 


=*  ?“=“■&[ T3T(?+TEiF+  (5^]  026  J 

|fci|  —  M** 


so  that  5=0, 

If  the  initial  condition  of  ^  and  the  potential  l/  are 
*  smooth  functions  of  x,y,a,  for  a  reasonable  time  interval 
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T  of  physical  interest,  then  the  local  truncation  error 
may  be  assumed  to  be  about  the  same  order  as  that  of  the 
round  off  error.  The  sufficient  condition  for  utility 
is  simply 

a*  &  _ _  028) 

Actually,  there  is  an  upper  limit  for  T  depending  on 
the  smoothness  of  V  as  well  as  on  the  boundary  and  on 
the  initial  condition  of  if  . 

(2).  Two  level  implicit  difference  scheme  (backward  time 
differences ) 

From  (123) 


+ 1  Pf  .  - 


+%i-e*P(i%  fc“v-'  ~*i  « > - 2  <TJ 

=  "  (129  ) 
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The  error  satisfies  the  equation 

An*l 


•  A.  e  r  H*l  I 

B«n~  T^a 

+ C)  4 &%P(-k  14  ^ )  ny^](~r} 

+  E£ *  C ■ Tfc  f?«„ b”'t ](4 

+  fd." 


=  £.  +  p" 

*  i A  j  ft 


:<jfl  . 


To  apply  (  117),  one  computes: 

RyLPi )  =  |  I  +  (  j£x>* +  (Sjjf  +  (tOtf)  | 

'I  £x>  +  £5o*  + 

I  -  iSr^’  (  tzi<?  4'  {&yf  ”*'  . 

For  si  s=  =  loo  9  the  sufficient  condition  for 
is 

. ,  ^  _ 222 _ 


}« 

(130) 


(131  ) 
utility 

(132) 
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(3).  A  simple  multilevel  difference  scheme  (Leap-Frog 
Scheme ) 

C  =  +  ihn-%*w 

+ C.4*e^  £  *,  j4-  v*)  -2 

+f^AeXP^^j 

+  %.l4  **?  (--k  {i'Apl-  & 

+  %ft;  e%Ph7C(-  *>**)  - 

(133  ) 

Since  $=2  for  three  levels,  C*  has  a  very  simple  form 
namely 


C 


ft  i) 


( 134) 


Therefore, 

««- ,+^-”  l 1 1 + 

(Sx>  'f'  feay/  +  ^a)*-) 

~  1  +  “**  +  d>)  .  (135) 

Again,  for  b/»  jfo-ioo ,  the  sufficient  condition  for 
utility  is 


***-  Trzn 


WL 


^(taxy  toy* 4'  i2t5l) 


C 136) 
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Let 


t  ^expc-t^x1’)  ,  £  =  £  expc-iw-x*) 


%  e*f>(+  ^X°)  ,  tf^  =ssr exp(+  isp. xp)  ' 


Then 


H-if-W  +  ^pc^HQi,  +gft)% 

(140a) 

T$  +  uA'%  +  +  #*’) 

+;A+ 4?^j  + 4?^  ° 

(140b) 

7r  *  4f^  +  e*P(-&j&’{(j?  +*!*’}£ 

+fi^'  +^f  A)"^&  +  0 

(140c) 

t& + 4?  A%  +  e*p(-  *J{I& + ifA') 

(l40d) 

(a).  Two  level  explicit  difference  scheme 

(forward  time 

differences ) 

<  r-ntl  _«*  * 

•  tWJ®  =  tip,  -  wivstp* ) 

—  It  P.K 

+(C£if  - €-;<)  + 

- > (  Klifrtrfyij* &+*£..<; r"  )  1  x. 

(141a) 
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Following-,  the  same  procedure  as  for  the  Schroedinger 
equation,  one  can  show  that  for  the  sufficient 

condition  (  HI  )  for  utility  is 


AX  ^ 


042  ) 


Naturally,  there  is  an  upper  limit  for  /  as  there  was 
for  the  Schrodinger  equation* 

(b).  Two  level  implicit  difference  scheme  (backward  time 
differences) 


(j  "W  .t  mi 

a r*(Af4)  / 

Tn#** 

-i(  Sk ia |t i- « + g/f <Ci «JK  =° C j « 


Ci«>+-sKC>*'+  exp**-*®*;  • 

\  2AX'  ^  J 

■ .  ■  /  tf"*  -  a:”*'  hH 

4*  /  (  +  J.e.A'$,  .  O 

2  4X*  -ftc” 


(143a) 


-mi  ^»*i*-i  H 

-(  &%„>]**■=  %.,*>  (143b) 
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■^c.  J  «,  -  fia>  +w  <$t  0  ) 

-  *  (  5ttt^.Saaffl  4  Ci  »irf°  €  i 


4) 


(143c) 


4  &%.>iX' +  «?«****“>  • 

f  /  Tnfl  .r«« 

-(  %^^Xh±0  +  Jetfqf^jAX’-  <%. .  w  (143,0 

Following  the  same  procedure  as  for  the  Schrodinger 
equation, one  gets 

i  .  i 


Therefore,  the  suff>  jnt  condition  for  utility  with 


(144  ) 


N~£,  =  loo  is 
ax 


(145) 
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(c)*  A  simple  multilevel  difference  scheme  (Leap-Frog 
scheme) 

•  +  iS/td*  .  'I 

\  2AXJ  KC^  ) 

+(%*vM>~}£±±Ltl  .f  i-S-A'i i"  „  ^ 

V  2dX'  frC  j  4)  / 

-*  (  -aii&4?.r  ±tlih*2  +  il A*  jr "  .  ))2AX*  ( 1 46a) 

2AX1  +Tc^  %<;j4)4r"A 

—*>♦1  r"*1 

**  ~  -  exp(^c^)  « 


fe'jo  +  i.i  v  ^ 

2  AX'  tfC  "  / 


^  _ 

r  h-l  m 

%;,»>  =  %it)  ~  ^cH’%..nAx'-  "MdUagGs i). 

-i( 

Cj«,  =  C;«  ’  TiA'C^x'-  . 


(146b) 


(146c) 


■) 


#  ,  JjLa'V*  • „ 

lV[  2 ax  +  tjc 

+a($Lu«*L1 $U&i)+Ae*(z''  ) 

*  2<ax‘ 

-  f  +  Al  a»ie*  )  ) 2/JX' 

2<jx*  Tc utsjM'J  OA  • 


(146d) 
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1  s  <; 


Then 

R  —  i  + 2  O&zjA'l  +  d47) 

=  i  ■*  Tax’ 

so  that  the  sufficient  condition  for  utility  with 
a f  =  =  ls 


3.  Won-Linear  Snlf-Focnslng  Equation 


Let  t  so  that 


(148) 


(149) 


(150a) 

(150b) 


(a).  Two  level  explicit  scheme  (forward  time  differences) 

■  Ci>  -  Cj>  - 

■+  tuioil  1  -f  — L  /  rAw^r) 

(^rj*  /ay  ^  2^r  ' 


4  L 

(^r;*  zf-ar  Hr  / 

+  CtA7/;)/+  . 


(151a) 


( i  51  b  ) 
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The  error, satisfies  the  approximate  relation: 


C = +-  StaLt&p!^ 


+  — 
a&Y 


+  2A C,,}** 4  +  *C-j) 


(152a) 


C<>“  Cjj  t{“  4-  E^>+fyf*EW 


+ 


-1~  -f  J.  a /aM  fir11 

\  2^»r  A  L(/LuiJ  +3(A,(;i))jEl(;i: 


+  u;Vj) 


(152b) 


This  yields 

*  -  I  1  +  zba*\  +  4  +  #f  +  *j«**' 

-  l  +  [W+  +  2  (3‘  -■  /  +  W-;  (153  ) 
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where 


i. - 5js I Aij, AiiJ .  j* 


For  a/2*  ~^tsslo°  ,  the  sufficient  condition  is 

for  all  1*1/1;  —  . 


(154) 


(b).  Two  level  implicit  scheme  (backward  time 
differences) 

■H  -W*l  HM  H«  .Ml 

A**1  AtHM  ~  Amm  , 

"r  (^r)fc 


242 

-r.M 


—  A 


l«J> 


A  +, 

„  A**  -A*". 

•'!(/>*•{)  ^  'W-ljf 

*4«j, 

n«;j)  + 

* 

(ar)‘ 

“  7a7 '( 

(O^ia. iCA<i"L)'Lr (A**1  ( k*"  4 

7CK 

II 

4>  • 

(155a) 


(155b) 

Following  the  same  procedure  as  in  the  first  two 
examples , 

*  <i«) 
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For  //*»  2^  =  /°°  ,  the  sufficient  condition  is 


a*  * 


(  157) 


(c).  A  simple  multilevel  difference  scheme  (Leap-Frog 
scheme ) 

a"*1  —  /f"  o 

nlMj>  nK.'^  *  1  ZdZ 

.  A*tA+lj}  +  . _ I  / hiU*!))  ~ 

‘  (d\r/  "r  2 &y  ) 

A"*  »  +  2i-  ALLi*0.~^JM 

+  ^«f)  i  jl  / ^i(iui) ~  Wl 

(ar)*-  mv  v  zar  / 

+pK,pU  (Cj, £}/&,•>}«*  .  ossb) 


(  158a) 


This  yields 


+2£^+r2-f  +  *(j, <}<)>*-  /  +  S<>*  .  (  159) 


<a  2  iter)1 


Therefore,  the  sufficient  condition  for  utility  of 
with  A /=*sr.z=(0o  is 


4/ 


(  160) 


for  all  "  =  /,£, - .  As/. 

Another  example,  the  laser-fluid  system  of  equations 
^considered  in  thin  study,  will  be  given  5n  the  next 
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section.  It  is  similar  to,  but  actually  simpler  than, 
the  self-focusing  example  just  considered. 
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PART  IV 

COMPUTER  SOLUTION  OF  THE  LASER-FLUID  EQUATIONS 

In  this  section  a  computer  solution  of  the  full  set  of  laser-fluid 
equations  is  presented.  The  two  major  difficulties  in  using  numerical  tech¬ 
niques  to  solve  differential  equations  by  computer  are  error  growth  and 
excessive  computation  time.  In  order  to  control  the  error  growth,  the 
utility  criterion  discussed  In  the  preceding  section  has  been  used.  Further¬ 
more,  highly  accurate  seven-point  difference  quotient  representations  of  the 
differential  operators  were  employed  to  reduce  truncation  error.  A  procedure 
for  obtaining  such  representations  Is  given  in  Appendix  II.  In  order  to 
handle  the  economic  problem  of  large  computation  time,  a  certain  amount  of 
efficiency  is  Introduced  by  minimizing  the  amount  of  core  storage  required 
of  the  computer.  This  was  accomplished  in  part  by  using  overlaying  techniques 
to  store  several  pieces  of  information  at  the  same  site  in  the  computer.  Thus 
information  Is  stored  only  as  long  as  it  is  needed  and  then  is  replaced  with 
current  material.  The  computation  time  was  also  minimized  by  making  use  of 
a  nonuniform  grid.  The  seven-point  difference  relations  allowed  a  relatively 
large  grid  size  without  undue  truncation  error  and  the  nonuniform  grid 
permitted  a  greater  grid  density  in  the  region  of  special  interest.  Thus  an 
accurate  solution  could  be  obtained  with  a  minimum  of  computation. 

The  laser-fluid  equations  were  solved  in  the  near  field  region  of  a 
laser  pulse,  initially  gaussian  in  both  r  and  z,  propagating  through  air  at 
1  atm  of  pressure  and  at  10°C.  A  cylindrical  geometry  was  used  and  cylindrical 
symmetry  (no  dependence  On  the  angle  <j>)  was  preserved  at  the  price  of  dropping 
the  gravity  term  in  the  Navier-Stokes  equation.  Having  cylindrical  symmetry 
amounts  to  a  considerable  simplification  in  the  problem,  so  that  the  inclusion 
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of  the  free  convection  effects  due  to  gravity  was  not  attempted  in  this 
analysis.  The  problem  described  above  amounts  to  a  mixed  Initial -boundary 
problem.  The  Initial  configurations  of  the  laser  beam  and  the  fluid  are 
specified  subject  to  certain  boundary  conditions  at  r  =  0  which  must  be 
satisfied  at  all  times.  Furthermore,  the  boundary  condition  at  z  =  0  Is  time- 
dependent,  because  the  tall  of  the  gausslan  must  be  fed  Into  the  spatial 
region.  For  the  numerical  solution,  a  spatial  mesh  of  grid  points  or  stations 
Is  'Jsed  to  represent  the  rz-plane.  At  a  given  Instant  In  time,  the  values  of 
the  various  dependent  variables  are  obtained  at  all  of  the  stations.  The 
difference  equations  are  then  employed  with  these  values  of  the  dependent 
variables  to  advance  a  step  In  time.  This  procedure  Is  repeated  over  and 
over  until  the  desired  time  Interval  has  been  traversed.  An  explicit  difference 
scheme  was  used  In  this  calculation,  because  such  schemes  are  simplest  to  handle. 

The  laser- fluid  equations  are  given  In  equations  (25a  )-(29 )  In  Part  I. 

As  mentioned  above,  the  gravity  term  was  dropped.  Also,  the  thermal  conduc¬ 
tivity,  k,  was  taken  to  be  constant  because  its  derivatives  are  very  small. 

The  equation  of  state  was  taken  to  be  the  Ideal  gas  law.  The  numerical  values 
used  for  the  various  parameters  are  the  same  as  those  given  for  the  linearized 
analysis  In  Part  II ,  because  the  same  temperature  and  pressure  were  used  for 

the  undisturbed  medium.  The  laser  frequency  w  and  the  dimensionless  absorption 

14  -1 

constant  aQ  were  chosen  to  be  1.773  x  10  sec  and  10,  respectively.  The 
wave  equation,  (25a),  for  linearly  polarized  light  In  an  absorptive  medium  Is 
taken  to  be 

2 

c2V2E  =  ^  (eE)  +  cxC^L  e).  (161) 

3t 

This  equation  is  an  approximate  equation  describing  an  electric  field  which  Is 
polarized  linearly.  Strictly  speaking,  of  course,  Maxwell's  equations  do  not 
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allow  cylindrlcally  symmetric,  linearly  polarized  beam. 


The  solution  of  (161  )  Is  taken  to  be  in  the  form 

i  iUt  -  k,z)  -  f  z 

E  -  ^  (E1  ♦  1E2)  e  ^  L  e  £  +  c.c.  ,  (162) 

where  E.|  and  E2  are  slowly  varying  functions  of  and  t  and  the  laser  frequency 

and  wave  number  are  related  by 

e0^>L  =  c  k^  ,  (l  62a ) 

where  eQe  Is  defined  in  (30f).  Substituting  (l 62  )  into  (161  )  and  dropping 
the  second  derivatives  of  E^  and  E2  with  respect  to  time,  one  can  put  the  wave 
equation  in  the  form 


f\*l+  c2  f„2E  +„^2^po  a2  1 

•i[rte-  +  s[r  [V  e2  taTt^S~'  /  te2| 

+  ^7  [“L  <^o  •  &  E1  -  ^!|lf  Ez] + 

+  &  (e  -  eo>  E2  ■  r  Ip  It  Ei + 


V  d?W  2 


.  J d2e  (zr\  c  _  L  de  /aA  ^2  _  d%\ 

2u^e  5?  \dt J  LZ  ”  o^e  l  Bp  \dtj  3z  ~  3z2  J 

-  s[r  w  -  sf?  pEi +  •  w-  -  ’) ♦  t  Ei] ♦ 

+  [“l  ^  e2  +  It  It  Ei]  + 

'  H  (e  '  Eo>  E1  '  llflt  E2  + 


1  \A  r  L  d£  3£  ^ 

■wLe  [d7  vV  E]  "  V  V  * 3t  32  "  WL  ai1"/! 


(163  ) 


(164  ) 
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Introducing*  now,  cylindrical  coordinates  and  using  cylindrical  symmetry 
and  the  other  simplifications  mentioned  above,  one  obtains  the  laser-fluid 
equations  In  the  following  form: 


(l  65  ) 


31 

3t 


1 1 1  f  1  3T  sh\ 

%  I  II?  7  »  1?) 


+  fsz 


<Ef  ♦  Ei> 


-az 
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0 


(168) 


3E1  _  c2  f  o  ..  3E1  .  3  E2  .  1  3E2  .  3  E2  .  Eo  -  e  3E2  .  .  1 

^ *7?" fo^V^)^-  +  TE2j 


.  ac 
+  1E 


"eo  - 

a;+ 


E-,  - 


j£  E 


>1  2e 


St  -2j  -  ^  <E  "  eo>  E2  "  Fe'  ftEl 


(170) 


£  3E,  2  I 

7^^EJ 


C 
U^E 


i/2e,  ia!!i  c2|cl  32eA1 


(171) 
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where  the  notations 


and  e"  h  & 
dp 


(172a) 


de 
=  dp 


have  been  employed  and 

e  =  eL  +  \  e2  (E2  +  l\)  (172b) 

The  easiest  way  to  obtain  a  variable  grid  size  Is  to  Introduce  a  trans¬ 
formation  to  a  new  Independent  variable.  Thus,  In  order  to  have  more  grid 
points  In  the  region  of  special  Interest,  small  r,  the  nonlinear  transformation 


(173  ) 


can  be  employed.  The  scale  value,  rQ,  can  be  chosen  later  according  to  the 
dictates  of  convenience.  After  this  transformation,  equations  (165)  -  (171) 
become,  for  x  j*  0: 


*k 


(l-x)2 
r_  3x 


M  v  +!!z 

rQx  r  3z 


+  -0~XJ  v  +  v 

r0  r  3x  vz 


(174) 
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(*|i  +  lh*ii„  !!i+v  !!z) 

(M  3z  rQ  V  ax  vz  ax  j 


(177) 
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(178) 
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(179) 


Because  of  the  symmetry  of  the  problem  and  the  regularity  of  the 
differential  equations,  the  following  boundary  conditions  must  be  satisfied 
at  x  =  0: 


3E2  _  3T  _  3p 

dX  3X  3X 


and 


( 180  ) 
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Equations  (174)-{179)  are  not  useful  at  x  ■  0  because  Indeterminate  ratios 
such  as  vr/x  appear  In  these  equations.  Applying  equations  (180  )  and  using 
L'HopItal's  rule  to  evaluate  the  Indeterminate  ratios,  one  obtains  the 
following  equations  for  use  at  x  ■  0: 
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Equations  ( 174)-( 186)  are  the  basic  equations  to  be  used  for  the  solution 
of  laser-fluid  problems  In  the  near  field  region  when  cylindrical  symmetry 
pertains.  These  equations  must,  of  course,  be  converted  to  difference  equations 
before  the  computer  solution  can  be  attempted.  The  difference  equations  will 
not  be  presented  here,  because  they  are  Included  In  the  computer  program  shown 
In  Appendix  C.  The  Interested  reader  will  be  able  to  locate  the  difference 
equations  In  the  program  listing.  As  Indicated  earlier,  seven-point  difference 
quotients  were  used  to  represent  the  differential  operators  appearing  In  (174)- 
(186).  These  difference  quotients  and  a  method  for  deriving  them  are  shown  In 
Appendix  B.  These  expressions  will  be  used  In  the  consideration  of  the 
utility  criterion  for  the  laser-fluid  system  of  equations. 
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Once  equations  (l 74)-(l 86)  have  been  obtained.  It  Is  vital  that  the 
utility  criterion  or  some  other  criterion  be  applied  to  determine  useful 
time  step  sizes  and  corresponding  grid  spaclngs.  In  spite  of  the  complexity 
of  these  equations,  the  utility  criterion  Is  extremely  easy  to  obtain.  This 
Is  one  of  the  attractive  features  of  the  utility  approach.  To  Illustrate  the 
ease  of  application  of  the  handy  sufficiency  tests  for  utility  described  In 
Part  III,  the  derivation  of  the  utility  criterion  for  the  laser-fluid  system 
will  now  be  given.  In  the  Interest  of  brevity  and  to  emphasize  the  simplicity 
of  the  derivation.  It  will  be  presumed  that  the  reader  Is  well  acquainted  with 
the  discussion  of  the  utility  method  given  In  Part  III  and  with  the  examples 
presented  there. 

An  examination  of  equations  (174)-(186)  reveals  that  It  Is  sufficient 
to  consider  either  (178)  or  (179)  and  Ignore  all  the  other  equations.  These 
two  equations  dominate  the  utility  criterion  and  either  (178)  or  (179)  can  be 
used,  because  either  choice  produces  the  same  condition.  Choosing,  then,  to 
deal  with  (178)  and  keeping  only  the  most  Important  terms,  one  can  obtain  a 
utility  criterion  from  the  equation 


3E1  c2 

aT***  25[e 


+ itsii  (,_&)  !fa 

rJT  ' s* 

o  o 


+ 


2k. 


if] 

l  dZ 


A  ^  2 
c  \  rt 


(187) 


2  2  2 

Using  the  fact  that  c  ■  eQe\  and  putting  e  *  eoe  «  1  one  can  simplify  (187) 
to  produce 


3E1  c 


(1-x)4  3^E2  .  (1-x)3  (l-2x)  3E2  3E1 

3x  x  3x  ^  3z 

L  o  o 


(188) 
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The  corresponding  seven-point  explicit  difference  equation  Is 


Fn+1 


*  E 


cAt 


»•  i  vuw 

l,jk  2k^ 


(189) 


where  the  difference  operators 


are  shown  In  equations  (B12) 


and  (B121)  of  Appendix  B,  respectively.  The  superscript  n  Is  the  time  Index 
and  j  and  k  are,  respectively,  the  x  and  z  Indices  for  the  spatial  grid  points. 
Thus, 


x  =  jAx  and  z  =  kAz 

The  various  coefficients  are  maximized  for  the  choices 
(l-*)-l  and 

where.  In  the  denominator,  x  =  jAx  -*■  Ax  because  (189)  does  not  apply  at  j  =  0. 
[Equation  ( 189)  was  taken  from  (178),  whereas  ( 185)  Is  the  appropriate  equation 
at  x  =  0.]  After  making  these  replacements,  one  obtains 


„  En 

■l.jk  M.jk 


CAt  j 

2\\ 


r02(Ax)2 


*j  + 


n  2kL  .  En  \ 
2,jk  “  Az  Ak  fcl  ,jk  j ' 


(190) 


where 


operating  on  subscript  j 

operating  on  subscript  j 


operating  on  subscript  k. 
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One  finds  that  the  precondition  (104)  Is  already  satisfied  because  no  diagonal 

jv 

term  arises  from  A,,  when  the  seven-point  scheme  (  rr  .  Is  used.  Thus  condition 
k  «  7,4 

(105)  can  be  used.  Using  equations  (B12)  and  (B121),  IT  Is  trivially  computed 
from  (190)  and  one  obtains 


R  =  1  +  6At, 


where 

I 


2k,_r02(Ax)2 


[  <|  +  f) 


/  3  3\ 

(7"  V 


h*  (i5* 


(191) 


then,  using  (116)  with  qQ  * 


kLro 


12 

T 


(Ax)‘ 


4,  one  obtains  the  utility  condition 

jl|. 

AZ I 


(192) 


This  Is  such  a  strong  constraint  that  one  Immediately  wonders  If,  perhaps, 
more  generous  sufficiency  theorems  for  utility  should  not  be  sought.  The 
utility  philosophy  suggests  a  more  practical  approach:  Try  the  condition  and 
then  try  to  violate  it  and  compare  the  results.  This  was  done  for  the  problem 
discussed  In  this  report  and  no  escape  from  (192)  was  possible.  In  fact,  If 
the  criterion  was  violated  by  a  factor  of  order  5  In  At,  then  classic 
Instability  phenomena  were  observed  in  the  computer  output.  Thus,  by  a  stroke 
of  bad  luck.  It  appears  that  (192)  must  be  obeyed. 


In  order  to  emphasize  the  Implications  of  (192)  for  the  study  of  the 
propagation  of  laser  pulses,  a  description  of  the  accessible  parameter  regime 
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will  now  be  given.  As  a  starting  point  for  this  discussion,  the  parameter 
values  used  in  the  actual  calculation  will  be  listed.  The  electric  field 


at  t  = 


where 


0  was  taken  to  be  of  the  form 

-4(f-)2  -4^o)  2 

E,  =  F  e  0  e  '  0 ' 


(193) 


e2  =  0  , 


full  ^  -width  of  E.|(r,z  ■  z  ,  t  **  0) 

•  [full  j  -width  of  IL  at  z  =  z  ,  t  =  0] 


IL  =  [on-axis  intensity  (in  ergs/cnr  sec)  at  r  =  0,  z  =  zpQ, 

t  =  0,  time  averaged  over  several  optical  periods] 

=  1  *£“  c  F2 
i  oe 

z  =  full  -  -width  of  E-|(r  =  0,z,  t  =  0) 


z  =  location  of  the  peak  at  t  =  0 


po 
F2  = 
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16  P, 
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— _  "  _  IS  -  ~ 

c  c  r  n 
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-  c  r 

oe  0 


eoe  ro  zo 


(194) 


=  peak  value  of  the  electric  field,  squared 

PL  -  total  power  (in  ergs/sec)  of  the  pulse  at  z  =  zp, 
t  -  0,  time  averaged  over  several  optical  periods 

U  =  total  energy  in  the  pulse  at  t  *  0,  time  averaged  over 
several  optical  periods. 


The  values  taken  for  these  quantities  were 
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r  =  200  cm 
0 

5 

zQ  s  9  x  10  cm  =  9  km 

a  13.5  x  10®  cm  =  13  km 
po 

F  =  4.7  x  10^  [erg/cm2]^2  (195) 

IL  =  3.3  x  10^  [ergs/(cm2  sec] 

PL  ■  5.2  x  1021  ergs/sec 
U  -  9.8  x  1016  ergs 

and  the  air  was  taken  to  be  Initially  In  Its  unperturbed  state  at  1  atm 
pressure  and  at  10°C. 

The  spatial  grid  was  composed  of  80  x  16  -  1280  mesh  points.  The  z- 
axls  was  evenly  divided  Into  80  steps  of  size  Az  *  0.45  x  105  cm  *  0.45  km 
beginning  at  z  =  0  and  extending  to  z  *  35.6  x  105  cm  *  35.6  km.  Thus  the 
peak  of  E-|  was  Initially  located  at  the  30th  mesh  point  on  the  z-axis  and 
its  1  -width  extended  from  the  20th  to  the  40th  mesh  point.  The  radial  variable 
x  has  the  range  0  <_ x  <_  1  and  this  range  was  evenly  divided  into  22  steps  of 
size  ax  -  |p  but  only  the  16  sites  closest  to  the  z-axis  were  used.  The 
more  distant  sites  correspond  to  radial  distances  greater  than  5  beam  half¬ 
widths.  The  first  step  away  from  the  z-axis  corresponds  to  the  radial 
distance  Ar  *  2-^yCl^  *  9.5  cm  s  (radial  half-width).  The  time  step  size 

was  taken  to  be  At  =  1 0” ?  sec  and  100  steps  were  made  so  that  the  time  Interval 
-5 

0  <  t  <  10  sec  was  trc'ersed. 

Taking  these  grid  sizes  and  time  steps  and  substituting  Into  the 
utility  condition  (192).  one  gets 
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(196) 


c"T  =  3M)  -  '40800  +  3536 


12  8_ 

^/(Ax)2  Az 


Thus,  for  the  chosen  step  size,  Az,  the  value  used  for  At  would  violate  the 
condition,  were  It  to  be  doubled.  Of  course,  the  condition  (192)  Is  only 
approximate,  but,  as  mentioned  above,  good  solutions  could  not  be  obtained 
for  At  -  10“6  for  Az.  Clearly  the  value  -gg  used  for  Ax  does  not  saturate  the 
Ax  piece  of  (196)  and  one  could  probably  use  Ax  as  small  as  Such  a  small 
step  size  for  Ax  would,  however,  require  three  times  as  much  spatial  mesh 
points  and  would  exceed  the  storage  capacity  of  the  computer  which  was  used. 


The  desire  Is  to  use  as  large  a  value  of  At  as  one  can.  In  this  regard, 
the  Ax  piece  of  (196)  Is  generous  and  would  permit  At  ~  10"6.  The  Az  step 
size  would  have  to  be  Increased  to  Az  ~  5  km  to  allow  this,  however.  Such  a 

large  step  size  would  be  larger  than  the  4.5  km  half-width  of  the  pulse 

selected,  so  that  no  details  of  deformation  of  the  pulse  could  be  observed. 

If  the  beam  Is  made  narrower  In  radial  extent,  the  Fresnel  length 

decreases  and  diffraction  effects  become  Important.  The  Fresnel  length  Is 

310  km  for  r  =  200  cm,  so  that  one  would  become  Involved  with  far  field 
o 

effects  If  the  beam  radius  were  decreased  by  more  than  a  factor  of  5. 

Making  the  pulse  longer  In  the  z-dlrectlon  expands  the  time  scale  over  which 

Interesting  effects  may  be  studied.  If,  on  the  other  hand,  the  pulse  Is 
shortened  In  the  z-dlrectlon,  then  one  must  shift  to  smaller  values  of  Az 
In  order  to  be  able  to  follow  details  of  the  development  of  the  pulse. 

Shifting  to  smaller  Az  requires,  because  of  (l 96 ) >  that  one  use  smaller 
values  of  At.  The  net  effect  Is  that  no  profit  Is  derived  from  using 
shorter  pulses,  because  they  can  be  followed  only  for  correspondingly  shorter 
times. 
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One  aspect  of  the  parameter  regime  has  not  yet  been  discussed*,  the 
range  of  power  for  the  beam.  Since  the  power,  P^,  does  not  appear  In  the 
utility  criterion,  Its  role  must  be  determined  by  experimentation  with  the 
computer  program.  Very  small  powers  are  not  Interesting  because  there  Is 
very  little  Interaction  with  the  fluid.  In  order  to  see  Instabilities  and  non¬ 
linear  effects  during  short  times,  one  would  wish  to  consider  beams  with  large 

power  densities.  The  extremely  large  values  shown  In  (195)  produce  Interesting 
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effects,  In  a  time  Interval  of  10  seconds.  Such  beams  can  not  be  followed  for 
more  than  about  100  time  steps,  however,  because  the  various  dependent  variables 
begin  to  develop  large  curvatures  and  vary  on  a  scale  smaller  than  the  meshT'' 
sizes.  Thus  If  one  wishes  to  follow  the  development  for  a  longer  period  of 
time,  the  mesh  sizes  must  be  decreased  and  eventually  the  time  step  will  have 
to  be  smaller,  and  then  many  more  time  steps  will  be  required.  In  this  regard, 
one  must  keep  In  mind  that  If  the  mesh  size  Is  decreased,  while  the  Initial 
pulse  size  Is  not  decreased,  then  more  mesh  points  will  be  required  and  the 
storage  capacity  of  the  computer  also  becomes  a  limiting  factor.  The  final 
remaining  option  Is  to  Increase  the  power  In  the  beam  even  more.  The  net 
effect  Is  that  the  large  curvatures  develop  faster  and  the  development  can 
be  followed  only  for  shorter  periods  of  time. 

One  final  remark  about  numerical  solution  of  the  laser-fluid  equations 
will  be  made  before  discussing  the  results  of  the  computer  calculation.  Strong 
growth,  Instabilities  and  nonlinear  effects,  cep  often  not  be  followed  because 
of  the  mesh  sizes  employed.  If  these  strong  oscillations  or  secular  growths 
are  generated  by  tiny  rapidly  changing  terms;  that  is,  If  the  instabilities 
arise  due  to  ripple  effects  which  become  strongly  enhanced,  then  a  crude 
mesh  size  can  smooth  these  effects  out  and,  thereby,  prohibit  the  occurrence 
of  the  strongly  growing  phenomena  by  removing  their  source.  Very  strong 
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Instabilities  were  found  In  Part  II  for  the  linearized  laser-fluid  equations. 
The  strongest  of  these  Instabilities  are  generated  by  very  short  wavelength 
ripple.  The  mesh  size  employed  In  the  present  calculation  will  begin  to 
wash  out  ripple  about  an  order  of  magnitude  larger  In  wavelength  than  the 
ripple  which  Is  most  strongly  amplified  In  the  linearized  analysis.  Thus, 
one  must  bear  In  mind  that  some  physical  sources  of  pulse  distortion  will 
be  excised  by  the  mesh  selected. 

Accepting  the  many  restrictions  noted  above,  we  have  examined  the 

17  -5 

propagation  of  a  200  cm  by  9  km  pulse  with  10  ergs  for  10  seconds.  The 
pulse  moves  three  kilometers  during  this  time  and  It  Is  possible  to  observe 
the  onset  of  the  laser-fluid  Interaction  In  some  detail. 

The  results  of  the  calculation  are  presented  In  Figures  [1]  -  [17]. 

The  electric  field  Is  conveniently  considered  In  terms  of  the  quantity 

|E|  s  [Ef  +  E*]  .  (197) 

where  E-j  and  E2  are  the  slowly  varying  electric  field  amplitudes  defined  In 
(162  ).  The  Instantaneous  electric  field  Is  thus  given  by 

E  =  | E  |  cos  [o^t  -  kLz  +  6e]  (198) 

where  the  phase  6E  Is  given  by 

6e  =  Tan”1  (E2/E] ) 

As  shown  In  (193).  at  t  =  0,  E2  Is  taken  to  be  zero  and,  consequently,  Is 
zero  Initially.  Thus 

|E|  =  E1  at  t  =  0  (199) 
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and  Ei  is  described  by  equations  ( 193)-( 195)  Initially.  This  Initial  pulse 
shape  Is  exhibited  In  Figure  [la]  and  In  Figure  [3]. 

In  Figure  [1]  the  z-profile  of  the  pulse  Is  shown  at  the  Initial  time, 

at  IQ’5  second,  and  at  two  Intermediate  times.  For  t  t  0,  the  pulses  are  not 

absolutely  symmetric  about  their  peaks.  In  order  to  exhibit  this  asymmetry, 

the  curves  are  plotted  as  a  function  of  ( z  -  z  |  ,  where  zQ  Is  the  center  of 

the  pulse.  This  device  allows  direct  comparison  of  the  leading  and  trailing 

edges  of  the  pulses.  The  center,  zc,  Is  defined  to  be  the  point  equidistant 

from  the  leading  and  trailing  edges  at  )E|  ■  1  (  V^2.  These  values  are: 

\  cnr  / 

t  *  0  t  *  6  x  10"®  sec  t  =  8  x  10"®  sec  t  =  10"5  sec 

(200) 

z_  *  30  z  ~  34  z  =  35.3  z„  ~  36.6 

c  c  c  c 

where  for  convenience,  distances  along  the  z-axis  will  be  given  In  units  of  the 

5 

grid  size:  Az  =  0.45  x  10  cm  =  0.45  Ion.  One  notices,  therefore,  from 
(200)  that  this  pulse  center  propagates  at  the  velocity  v£  *  2.97  x  1010  cm/sec, 
the  velocity  of  light.  The  pulse  peaks,  however,  are  observed  to  drift 
backward  with  respect  to  zQi 

t  =  0  t  =  6  x  1  O’®  sec  t  =  8  x  1  O’®  sec  t  =  10” 5  sec 

„  (201) 
z  *  30  z  s  34  z  s  35.2  z„  z  35.3 

P  P  P  P 

so  that  after  10’®  seconds,  the  peak  has  lost  about  two  thirds  of  a  kilometer 
with  respect  to  zc>  Note  that  the  exponential  damping  factor  shown  In  (i 62  ) 

Is  not  Included  In  the  quantity  |E|  appearing  In  the  graphs.  For  air,  this 
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factor  Is  larger  than  0.95,  even  at  t  ■  10  sec.  Other  than  this  effect,  very 
little  energy  Is  lost  from  the  beam  due  to  heating  of  the  fluid,  so  the  dis¬ 
tortion  effects  shown  In  Figure  [1]  are  rather  minor  and  are  noticeable  only 
near  the  peak  of  the  pulse.  Extra  detail  of  this  peak  distortion  Is  shown  In 
Figure  [2].  For  purposes  of  this  display,  the  leading  edges  have  been  placed 

-93- 


3/  \l/2 

together  so  that  the  curves  Intersect  at  |E|  =  10°  and  the  corres- 

Vcmv 

ponding  abscissa  has  been  labeled  42.5,  the  location  of  this  point  at  t  =  0. 

The  retrograde  peak  motion  and  the  corresponding  loss  of  fore  and  aft  symmetry 
In  the  vicinity  of  the  peak  are  plainly  seen. 

The  radial  beam  profile  Is  exhibited  In  Figures  [3]  and  [4].  The 
radial  slice  shown  at  each  value  of  the  time  Is  taken  through  the  position  of 
the  maximum,  zp,  In  the  z-proflle.  In  Figure  [3]  the  entire  beam  profile 
Is  shown  for  the  Initial  and  final  times  only.  Comparison  of  the  curves 
reveals  a  small  on-axis  Increase  extending  out  to  the  beam  halfwidth  (half  of 
the  full  £  -width)  at  -j-  *  100  cm.  The  effect  amounts  to  a  47%  Increase  In 
the  on-axis  Intensity.  Details  and  Intermediate  states  are  given  In  Figure  [4]. 

Although  the  radial  peak  Is  on  axis  In  the  slice  through  the  peak  In 

the  z-proflle,  this  Is  not  the  case  for  slices  taken  behind  zp.  For  example, 

at  t  *  10”5  sec,  the  principal  peak  Is  at  zp  =  35.3.  As  one  moves  away  from 

this  peak  toward  the  trailing  edge,  the  radial  peak  moves  off  axis  giving  a 

maximal  effect  near  z  *  31.  The  principal  peak  of  the  pulse  Is,  however, 

always  on  axis.  The  radial  profile  at  z  ■  31  Is  sh*wn  In  Figure  [5].  This 
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Is  clearly  only  a  small  detail  at  t  ■  10  sec.  The  position  of  this  off- 
axis  secondary  peak  Is  also  located  In  Figure  [9]  and  marked  with  tiny  squares. 

In  order  to  follow  the  development  of  the  phase  of  the  electric  field, 
the  quantity  |Ei |  Is  plotted  In  Figures  [6],  [7],  and  [8],  for  t  f  0.  Of 
course,  at  t  »  0,  ]Ej|  *  |E|  and  the  phase,  5^,  Is  zero.  In  these  three  figures, 
the  graph  of  |Ej  Is  marked  with  dotted  lines  for  comparison.  The  corresponding 
value  of  | E2{  can  be  deduced  from  these  figures,  using  (197).  These  figures 
show  far  more  dramatic  effects  than  the  curves  discussed  above.  At  t  * 

6  x  10"6  sec,  the  phase  Is  still  nearly  zero  and  E1  Is  positive  everywhere. 
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At  t  *  8  x  10"®  sec,  however,  has  changed  sign  over  a  three  kilometer 
region  extending  from  slightly  In  front  of  the  peak  of  |E|  toward  the  trailing 
edge  of  the  beam.  This  Is  clear-cut  evidence  of  the  onset  of  laser-fluid 
Interaction  In  the  trailing  edge  of  the  beam.  It  Is  clear  that  the  front  of 
the  pulse  and  the  distant  tall  are,  as  yet,  unaffected  by  this  Interaction. 

It  Is  Interesting  that  the  unperturbed  part  of  the  leading  edge  does  not 
reach  as  far  back  as  the  principal  peak.  Thus,  the  peak  already  feels  the 
effects  of  the  Interaction  to  some  degree  and  comparison  of  (200)  and  (201) 
reveals  that  the  peak  will  now  begin  to  lose  ground  with  respect  to  the  center 
of  the  pulse.  This  effect  has  already  been  noted  In  the  graphs  of  |E|  . 

Since  there  are  now  places  where  |E^|  Is  zero.  It  Is  clear  that  the 
phase  goes  to  j  at  these  sites.  The  amplitude  E2  responds  strongly  at  those 
places  where  E1  *  0,  fulfilling  the  obligation  to  conserve  power.  One  notices 
that  the  graph  of  |E|  remains  very  smooth,  giving  no  indication  that  the 
phase  Is  varying  rapidly.  Figure  [8]  shows  the  later  development  of  the 
region  In  which  E^  changed  sign.  The  region  In  which  E-j  <  0  Is  now  seven 
kilometers  long,  nearly  as  big  as  the  ^  width  of  |E|.  This  region  has 
advanced  no»  to  a  point  well  In  advance  of  the  principal  peak  and  extends  back 
far  Into  the  tall.  It  appears  that  this  node  is  propagating  forward  at  nearly 
four  times  the  speed  of  light.  Furthermore,  there  has  been  another  sign 
reversal  of  E^  slightly  behind  the  peak.  It  Is  this  kind  of  oscillatory 
behavior  In  E^ ,  with  large  variations  on  the  scale  of  the  chosen  mesh  size, 
that  brings  a  halt  to  further  observation  of  the  beam  development  by  this 
method. 


On  the  rz-plane  shown  In  Figure  [9],  the  constant  phase  curve,  E^  =  0, 
Is  shown  In  detail  at  t  ■  10"S  sec.  Also  marked,  with  small  open  circles.  Is 
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the  on-axis  extent  of  the  similar  curve  at  t  *  8  x  10“®  sec,  encountered  In 

Figure  [7].  Also  Indicated  on  the  same  figure  Is  the  locus  of  off-axis 
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maxima  In  the  z-proflle  of  the  pulse,  both  at  t  ■  0  and  t  a  10  sec.  At 

t  ■  0,  the  pulse  Is  described  by  (193)  and  clearly  the  off-axis  maxima  In 

z  are  all  positioned  at  zpQ  =  30  .  As  the  pulse  propagates,  however,  the 

peak  moves  slower  than  the  velocity  of  light,  as  previously  noted.  The  off- 

axis  portions  of  the  pulse,  however,  have  much  smaller  Intensity  and, 

consequently,  interact  very  little  with  the  fluid.  These  portions  of  the 

pulse  will  suffer  no  delay,  and  move  steadily  ahead  of  the  principal  peak. 

One  notices  that  at  two  thirds  of  the  radial  halfwidth  of  the  beam,  the  delay 

disappears  almost  completely.  As  mentioned  above,  the  small  squares  locate 
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the  secondary  peaks  present  at  t  *  10  sec. 


The  phase  front  information  presented  in  Figures  [6]-[9]  is  an  interesting 
feature  of  the  results  of  the  computer  solution.  Since  lE^  turns  out  to  be  a 
rapidly  varying  function  of  time,  It  Is  of  Interest  to  attempt  to  understand 
the  mechanism  responsible  for  the  behavior  of  E-j.  In  order  to  understand  this 
behavior,  one  must  realize  that  the  phase  depends  on  the  state  of  the  fluid. 

The  state  of  the  fluid  given  by  the  computer  calculation  Is  shown  In  Figures 
[10]-[16].  Before  these  figures  are  discussed,  however.  It  Is  convenient  to 
examine  certain  analytic  estimates  for  the  fluid  variables.  Such  estimates 

will  allow  Insight  Into  the  behavior  of  E^,  and,  later,  will  facilitate  the 

« 

discussion  of  the  computer  results  for  the  fluid  variables. 


In  order  to  describe  the  behavior  of  E^,  It  Is  useful  to  write  the 
electric  field  In  the  form 


/ew, 

1(w,  t - —  z) 

E  ■  £  (r,z-ct)  e  L  +  c.c. 


(202) 
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where 


„  9  4/Z-ct\2 

-4(f-)2 

*  _  F  .  ro  a  0 
€o  7  7  e  e 


e  -I 


( 202a) 


( 202b) 


where  F  Is  a  slowly  varying  amplitude,  Is  the  local  density  excess, 

(p-P0)»  and  zpo  has  been  put  equal  to  zero  for  the  present  considerations. 
Combining  equations  (202)-(202b),  the  electric  field  can  be  put  In  the  form  shown 
In  (162  )  with 


These  expressions  agree  with  (193)  at  t  =  0  and  offer  a  way  to  estimate  the 
behavior  of  E-j  at  subsequent  times.  On  the  basis  of  (203)  the  nodes  of  E-j 
might  be  expected  to  be  determined  by 
e  -1 

cos  (kLz  Ip  pl^  s  0  (204) 

Actually,  this  expression  should  be  modified  slightly  If  one  wishes  to  attempt 
to  get  quantitative  agreement  with  the  computer  solution.  It  is  clear  that 
(203)  requires  E2  to  vanish  at  z  =  0  at  all  times.  This  Is  not  the  same 
boundary  condition  which  was  used  In  the  computer  solution.  Actually  one 
should  use 
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E-|  ~  cos  \p 


wl  th 


dz'  . 


The  Integrand  In  (205)  should  be  evaluated  at  z1  and  t1,  where 


c(t'  -  t)  =  z‘  -  z. 


(205) 


Although  (205) will  be  used  In  the  future  to  analytically  describe  a  laser  beam 
quantitatively,  for  the  present  (204) may  be  employed  to  gain  Insight  Into  the 
behavior  of  E-|. 


It  Is  clear  from  (203)  and  (205)  that  the  mechanism  responsible  for  the 
behavior  of  E^  Is  easily  exhibited.  To  actually  follow  the  behavior  of  E.| , 
however.  It  Is  clearly  necessary  to  determine  the  state  of  the  fluid.  In 
particular,  the  density  excess,  p1 ,  must  be  obtained  as  a  function  of  time  and 
position.  In  order  to  analytically  describe  the  fluid  for  the  time  Interval 
and  parameter  ranges  of  the  computer  solution,  the  laser-fluid  equations  may 
be  simplified  to 


P 


o  at 


Y-l  3pl  .  ac  72 

»  w  ^ 


(206) 


and 


where  T1  Is  the  local  temperature  excess,  T-TQ.  One  must  recall  that  the 
o 

Intensity,  cE  ,  appearing  In  (206) Is  a  function  of  time  and  position. 
Integrating  (206)  from  zero  to  t  and  combining  the  resulting  equation  with 
(207) to  eliminate  T1 ,  one  can  obtain 


(207) 
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(208) 


It  Is  straightforward  to  Integrate  this  expression  and,  although  the  details 


will  not  be  given  here.  It  Is  clear  that  p.,  will  have  the  form 

?  1 
8r* 


F(z,t) 


(209) 


so  that  p  <  p  on  axis  and  there  Is  an  off-^axls  maximum  In  the  density.  In 
other  words,  there  will  be  a  pile-up  of  the  fluid  at  a  distance  r  ~  -j- 
from  the  axis. 


The  temperature  distribution  can  easily  be  obtained  by  Integrating  (206) 
from  zero  to  t: 


-8 


T.  ~ 


VPL 

cpocv 


r2 

7 


f 


-8(1-)' 

o 


z-ct 


d  v 
z_ 


(210) 


A  negligible  term  Involving  p<|  can  be  evaluated  using  (209)  and  has  been 
dropped  to  get  (210). 


Combining  (210)  with  (204)  the  behavior  of  can  be  visualized  and 
studies  analytically.  Closed  contours  such  as  those  shown  In  Figure  [9]  are 
predicted  and  other  qualitative  features  are  correct.  A  detailed  comparison 
of  this  analytical  procedure  and  the  computer  result  Is  In  progress  and  It  is 
now  clear  that  striking  quantitative  agreement  Is  obtained.  This  success  Is  of 
great  Interest  because  the  analytic  procedure,  unlike  the  present  computer 
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solution,  is  not  limited  to  10 


-5 


seconds. 
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The  computer  result  for  the  state  of  the  fluid  at  t  =  10  sec  will 
now  be  considered.  For  this  discussion,  It  is  useful  to  keep  several 
characteristic  distances  In  mind.  Since  the  state  of  the  fluid  Is  governed  by 
the  Intensity  profile 


IL(r,z,t)  =  \  ^  c  |E|2  , 


(211) 


rather  than  by  |F|  directly,  the  following  Initial  parameters  are  relevant. 
The  full  j  -width  of  the  z-proflle 


of  P,  =  —  =  14.1  In  units  of  Az  s  0.45  km 
*  2 


At 


if'zpi5. 


=  0  so  that 


3Z 


il 

9Z 


(212) 

(213) 


has  extrema  separated  by  10  units.  The  full  ^  -width  of  the  radial  profile 


cm 


of  I.  =  —  =  0.707  r„  *  141 

L  Jt  0  , 

r  3  i 

At  r  =  r  -  50  cm,  — J=-  «  0  , 

3r 


(214) 

(215) 


IL 

so  that  gp-  has  a  maximum.  The  maxima  on  opposite  sides  of  the  axis  are 
separated  by  100  cm. 


The  local  temperature  excess,  T  -  TQ,  where  TQ  =  10°C  Is  shown  In 

Figure  [10]  as  a  function  of  z  at  t  =  10’^  sec.  This  curve  Is  In  complete 

quantitative  agreement  with  the  analytical  result  shown  in  (210).  The  hottest 

place  in  the  beam  lies  on  the  axis  at  z  ~  33.3.  Since  the  intensity  peak  is 

at  zn  z  35.3,  it  is  clear  that  the  thermal  peak  is  lagging  behind  the  intensity 
P 

peak.  Since  z  ~  36.6  at  this  time,  it  is  clear  that  the  thermal  peak  is 
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almost  exactly  midway  between  the  Initial  and  final  pulse  centers.  Thus 
one  finds,  as  expected,  that  the  thermal  peak  propagates  at  velocity  |  for 
small  times.  This  and  other  properties  of  the  thermal  profile  are  readily 
understood  on  the  basis  of  the  following  considerations.  The  temperature 
responds  to  the  heat  deposited  In  the  medium,  so  that  from  (210): 


t‘z>  -  to = 

p  po' 


p“t;IL  -  0» 


(216) 


Due  to  the  symmetry  of  IL  and  the  fact  that  It  Is  gausslan.  It  follows  from 

(216)  that  T(z)  should  reach  a  maximum  midway  between  and  zp  for  small 

times.  Furthermore,  the  graph  of  T(z)  -  TQ  should  be  symmetric  about  Its 

maximum.  Both  of  these  features  are  evident  In  Figure  [10].  Since  the  peak 
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moves  much  less  than  Its  halfwidth  In  10  seconds,  the  Integrand  In  (216) 

Is  essentially  constant.  Evaluating  IL(r,z',0)  at  the  midpoint  of  the 
Interval ,  one  gets 


T(z) 


a 


I L(r,z 


0) 


(217) 


From  (217)  one  concludes  that  the  width  of  the  thermal  distribution  should 
equal  the  width  of  the  Intensity  distribution.  Indeed,  one  sees  in  Figure  [10] 
that  the  temperature  distribution  has  width  14.1,  which  Is  to  be  compared 
with  (212). 


The  corresponding  radial  temperature  distribution  at  z  *  33,  the 
position  of  the  maximum  at  t  ■  10"®  se.\,  Is  given  in  Figure  [II].  This  curve, 
also,  Is  In  complete  quantitative  agreement  with  the  analytical  result  shown 
In  (210).  The  temperature  has  reached  a  maximum  of  more  than  1000°K  on  axis 
and  the  full  width  of  the  distribution  is  found  to  be  0.67  rQ  *  134  cm,  5% 
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narrower  than  the  Initial  radial  Intensity  width.  One  might  have  expected 

the  temperature  distribution  to  be  broader  than  the  Intensity  profile 

because  the  large  radial  velocity  of  the  fluid  should  carry  some  of  the 

deposited  energy  away.  In  fact,  this  effect  may  possibly  be  observed  in 

the  following  way.  One  might  compare  the  temperature  distribution  not  to  the 

original  gaussian  Intensity  profile,  but  rather  to  the  ^  width  of  the  radial 
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intensity  profile  at  t  *  10  sec.  This  final  intensity  profile  has  a  width 
of  128  cm  or  5%  less  than  the  temperature  width.  The  average  of  the  two 
Intensity  widths  Is  134.5  cm,  almost  exactly  the  observed  temperature  width. 
This  average  may  be  the  best  measure,  because  the  thermal  peak  Is  midway 
between  the  initial  and  final  Intensity  peaks. 


The  z-component,  vz,  of  the  fluid  velocity  Is  shown  In  Figure  [12]  as 

a  function  of  z.  A  double  log  plot  Is  used  which  omits  values  of  vz  between 
-3  -3 

10  cm/sec  and  <*10  cm/sec.  This  kind  of  plot  allows  negative  values  of 
vz  to  be  plotted  below  the  "axis".  The  zero  of  the  velocity  distribution 
occurs  around  z  =  32.2,  so  the  "center  of  velocity"  lags  slightly  behind  the 
thermal  maximum.  From  (169)  one  might  expect  to  find 


v 


II  9Il 

Z  ~  3z  ~  3z 


(218) 


so  that  the  peaks  In  Figure  [12]  would  be  separated  by  10  units  according  to 
(213).  Indeed,  the  peaks  are  found  to  be  separated  by  10.2  units.  Furthermore, 
since  the  temperature  curve  Is  symmetric  about  Its  maximum,  (218)  would  suggest 
that  vz  should  be  antisymmetric  about  Its  zero.  This  effect  Is  certainly 
observed  In  Figure  [12].  The  velocity  distribution  Is  delayed  with  respect 
to  the  temperature  distribution,  but  this  symmetry  property  is  unaffected. 
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At  the  same  value  of  z,  corresponding  to  the  center  of  velocity  of  the 
z  component,  the  radial  component,  vr.  Is  plotted  In  Figure  [13].  This  value 
of  z  corresponds  also  the  the  largest  radial  velocities,  so  that  z  =  32  might 
also  be  termed  the  site  of  greatest  kinetic  energy  in  the  fluid.  The  radial 
velocity  maximum  Is  forty-six  hundred  times  larger  than  the  maximum  axial 
velocity.  This  effect  arises  because  of  the  great  disparity  In  the  Intensity 
widths  In  the  two  directions.  As  a  matter  of  fact 


(v  )  max 

—I - =  4625 

(vz)  max 


and 


9  x  105 

loo 


=  4500. 


The  curve  of  vr  Is  forced  to  go  to  zero,  as  r  goes  to  zero,  by  the  boundary 

conditions  shown  In  (180  ).  One  notices,  however,  that  the  peak  Is  located 

at  the  distance  0.25  rQ  from  the  axis,  exactly  the  location  of  the  maximum  of 
3L 

g —  shown  In  (215).  Thus  one  finds 


v 


3r 


II 

3r 


(219) 


as  would  be  expected  from  (177) 


Bringing  up  the  rear  In  the  sequence  of  effects  Is  the  density  minimum 
at  z  «  31.5.  The  density  decrement  -(p-pQ)  Is  shown  In  Figure  [14],  Using 
(207),  one  might  expect 


(220) 


so  that  p  and  T  will  have  the  same  z  dependence. +  There  should,  however,  be  a 
double  time  delay,  since  two  time  Integrations  are  Indicated  In  (220).  The 


^One  gets  the  same  conclusion  by  analyzing  the  function  F(z,t)  shown  In 
(209)  in  a  manner  analogous  to  the  reasoning  employed  to  obtain  (217). 
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width  of  the  density  decrement  Is  found  to  be  14. 5,  about  3%  wider  than  the 

Vp 

thermal  and  Intensity  widths.  The  maximum  fractional  decrement  ~ —  ~ 

-5  P° 

4  x  10 


The  radial  density  distribution  is  exhibited  In  Figure  [15]  at  z  =  32, 

the  location  of  the  density  minimum  In  the  z-proflle.  Again,  a  double  log 

plot  is  given  so  that  both  positive  and  negative  density  excesses  can  be 

conveniently  represented.  This  time  a  density  pile-up  Is  observed  because 

the  fluid  has  been  blown  away  from  the  axis  so  fast  that  a  compression  wave 

is  generated.  The  zero  In  the  graph  Is  at  0.36  rQ,  right  at  the  halfwidths 

of  the  thermal  and  Intensity  distributions.  Thus  Inside  the  thermal  halfwidth 

the  density  Is  depressed;  outside  the  fluid  has  piled  up.  The  radial  density 

profile  shown  In  Figure  [15]  is  In  excellent  agreement  with  that  predicted 

In  (  209).  For  example,  the  zero  observed  at  =  0.36  Is  predicted  to  occur 

ro 

at 


£-•  —  -  0.354 
ro  JB 


Similarly,  the  location  of  the  peak  observed  at 
at 


0.51  Is  predicted  to  occur 
ro 


Similarly,  the  ratio  of  peak  height  to  valley  depth  Is  also  correctly  predicted. 
As  a  matter  of  fact,  when  one  takes  the  trouble  to  evaluate  the  function  F(z,t) 
appearing  In  (209),  he  finds  precise  agreement  between  (209)  and  Figure  [15]. 
Thus,  both  (209)  and  (210)  are  in  complete  quantitative  agreement  with  the 
result  of  the  computer  calculation. 


Figures  [16]  and  [17]  exhibit  the  parade  of  effects,  illustrating 
graphically  the  various  delays,  pulse  shapes,  and  widths.  Physically  the 
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delays  make  sense.  First  the  beam  blasts  through,  heating  the  fluid  as  It 
passes.  As  explained  above,  the  temperature  maximum  moves  at  j  and,  thus, 
behind  the  laser  peak.  As  this  temperature  wave  passes  along,  the  fluid  picks 
up  kinetic  energy  and  the  flow  velocities  Increase.  The  center  of  this  effect 
trails  the  heat  wave,  allowing  time  for  the  fluid  to  respond.  Then,  as  the 
fluid  begins  to  flow  away  from  the  propagating  center  of  velocity,  density 
deficits  are  left  In  the  wake  and  corresponding  radial  compression  waves  set 
out  from  the  beam  axis. 

To  summarize  the  results  briefly.  In  portions  of  the  pulse  where  the 
Intensity  Is  small  there  Is  very  little  Interaction  with  the  fluid  and  these 
portions  move  without  appreciable  distortion.  The  peak  of  the  pulse,  however. 
Interacts  fairly  strongly  with  the  fluid  and  the  peak  Is  delayed  relative  to 
the  center  of  the  pulse.  A  parade  of  effects. ensues;  the  center  and  edges  of 
the  pulse  are  followed  by  the  peak,  which,  In  turn.  Is  followed  consecutively 
by  the  thermal  wave,  the  center  of  velocity*  and  the  density  waves.  The  front 
edge  of  the  pulse  propagates  without  appreciable  distortion,  but  strong  phase 
oscillations  are  set  up  near  the  peak  and  rapidly  overtake  the  undistorted 
front  section  Indicating  that  soon  the  entire  beam  will  be  distorted  to  some 
degree.  The  strongest  Instabilities,  predicted  In  the  linearized  analysis  of 
Part  II  did  not  appear  because  they  are  generated  by  ripple  with  wavelength  an 
order  of  magnitude  smaller  than  the  mesh  size  used. 

There  Is  very  little  hope  of  obtaining  computer  solutions  of  the  laser- 
fluid  equations  except  In  the  tightly  limited  regime  reported  here,  unless  a 
different  calculatlonal  procedure  can  be  devised.  Since,  however.  It  appears 
that  a  certain  amount  of  analytical  headway  has  been  made,  there  Is  reason 
to  believe  that,  with  appropriate  combination  of  analytical  and  computer  methods. 
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the  beam  can  be  followed  for  considerably  longer  periods  of  time, 
currently  being  directed  toward  this  objective. 


Effort  is 
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PART  V 


BEAM  FED  TURBULENCE 


Many  experiments  have  been  reported  for  which  theories  assuming  steady 
state  beam  profiles,  after  Initial  transients  die  out,  provide  rather  good 
explanations  of  the  principal  features.  However,  that  Is  probably  true  only 
because  these  experiments  are  conducted  at  relatively  low  power  fluxes. 
Theoretically,  one  expects  a  time-dependent  state  of  the  system  because  of 
the  Instabilities  discussed  earlier.  Such  Instabilities  are  not  observed  In 
practice  because  they  cannot  develop  within  the  distances  allowed  for  the 
propagation  of  beams.  However,  alternate  considerations  for  a  beam  of  finite 
cross  section  suggest  that  the  beam  may  drive  the  fluid  Into  a  time-dependent, 
or  turbulent,  state  at  powers  which  are  not  completely  unreasonable. 

It  may  be  Impossible  to  prove  analytically  that  such  a  turbulent  state 

develops,  because  the  Investigation  of  hydrodynamic  stability  Is  very 

difficult  even  for  the  simplest  flows.  However,  an  argument  can  be  made  from 

dimensional  considerations,  an  approach  that  promises  to  be  very  useful. 

Namely,  for  the  problem  of  a  beam  of  radius  a  and  power  flux  I  passing  through 

air.  It  Is  possible  to  estimate  a  parameter  W,  which  plays  the  role  of  an 

effective  Reynolds  number  for  our  problem.  It  will  be  shown  that  the  parameter 

W  takes  on  values  of  the  order  of  30,000  for  a  beam  with  Intensity  I  =  1  kilowatt/ 
2 

cm  of  radius  1  meter.  Since  It  Is  known  that  flows  with  Reynolds  numbers 
substantially  lower  than  30,000  are  turbulent,  the  flows  for  the  laser-heated 
atmospheric  path  should  also  be  expected  to  show  significant  time  dependence. 

Consider  the  equations  of  motion  for  the  air,  and  the  equation  governing 
heat  transfer,  which  take  the  following  form  If  It  Is  assumed  that  the  air  can 
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2  2 

be  assumed  Incompressible.  (That  amounts  to  dropping  terms  of  order  u  /v*, 
where  u  Is  a  typical  flow  speed,  and  v  Is  the  speed  of  sound.  For  the 
problems  under  consideration,  u  /v‘  will  be  less  than  10  ,  and  the  Incompres- 
slble  fluid  approximation  will  be  quite  good.) 

P  (v  •  7)  V  =  -v  P  -  8T,pg  +  n  v2  v  +  (n+n‘)  7(7  •  v) 

po  po 

p  (7  •  v)  =-(v  •  7)  p  ■  ep  (v  •  7)  T1  (221) 

P  cp  (v  •  7)  Ti  =  ^  +  7  •  (k7T] )  +  jj-  [n  (v^j  +  Vj^^  +  n'  (7  •  v)2] 

Assuming  the  Reynolds  number  Is  high,  the  Inertial  terms  will  dominate  the 
viscous  terms  In  the  Navler-Stokes  equation.  Thus,  there  must  be  a  balance 

between  the  Inertial  terms  and  the  bouyancy  forces,  which  Implies  that 

2 

pu  /a  ~  p3T-jg,  where  p  Is  the  density  of  air,  3  Is  the  coefficient  of  thermal 
expansion,  Tj  Is  a  typical  value  for  the  temperature  rise,  and  g  Is  the 

2 

acceleration  due  to  gravity.  The  pressure  variation  will  be  of  order  pu  . 

In  the  heat  transfer  equation,  the  convection  term  will  dominate  the  conduction 

term,  and  the  beam  heating  will  overwhelm  the  viscous  dissipation,  so  that  there 

must  be  a  balance  between  heat  deposition  from  the  beam,  and  convective  heat 

transfer,  which  Implies  that  pCp  uTj/a  **  ctl.  Combining  these  two  relations 
3  2 

we  find  that  u  **  a3a  Ig/pCp  .  With  that  expression  for  u,  we  then  define  a 
parameter  W,  which  Is  expected  to  Indicate  regimes  where  steady  flow  and  where 
time-dependent  flow  may  be  anticipated.  W  Is  an  estimate  of  the  relative 
Importance  of  Inertial  terms  to  viscous  terms  In  controlling  the  flow: 

aupn  on  ?  1/3 

W  3  3  a  ~  (<*3a  ig/pCp)  (222) 

10  2 

For  a  beam  with  I  3  10  ergs/cm  •  sec,  a  =  100  cm,  TQ  =  273°K  and 
a  =  3  x  10”6  cm"1,  W  =  30,000. 
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For  many  experiments  described  In  the  literature,  the  values  of  W  are 
much  smaller,  and  thus  one  would  not  expect  aqy  turbulent  fluid  flow  to  be 
observed.  For  example.  In  the  original  experiment  of  Gordon,  Lelte,  Moore, 
Porto,  and  Whlnnery^  the  parameter  W  takes  on  a  value  about  10"^,  and  In 
the  more  recent  experiments  of  Smith  and  GebhardtP^  W  Is  of  order  10. 

The  parameter  W  Introduced  here  Is  different  from  the  Grasshof  number, 
which  Is  referred  to  In  some  discussions  of  the  convective  flows  set  up  by 

n  q 

the  absorption  of  energy  from  a  laser  bearn.1^  In  fact,  the  conceptual  basis 
for  using  the  Grasshof  number  In  a  discussion  attempting  to  explain  the  tran¬ 
sition  between  smooth  flow  and  time-dependent  flow  seems  less  relevant  because 
the  Grasshof  number  appears  to  be  more  sensible  when  the  thermal  bouyant 
forces  are  balanced  by  viscous  forces.  In  the  present  discussions,  the 
thermal  bouyant  forces  are  balanced  by  Inertial  effects. 

We  are  planning  experiments  to  determine  the  critical  value  of  M,  VJ  , 

Ci 

which  determines  the  onset  of  turbulent  convective  flows  for  the  geometry 
appropriate  to  laser  beam  transmission.  It  Is  also  our  aim  to  attempt  a 
theoretical  evaluation  of  this  critical  value.  At  the  present  time  we  can 
only  speculate  that  Wcr  may  be  between  10^  and  10*. The  theoretical  approach 
appears  fairly  difficult  because  the  question  of  the  stability  of  flows  even 
without  heat  sources  has  only  been  answered  theoretically  for  very  simple 
geometries. [20-21 ]  yhe  qUest^on  0f  stability  for  fluids  which  are  heated  or 
cooled  appears  to  have  been  treated  mainly  for  cases  In  which  the  fluid  would 
be  motionless,  and  has  not  been  explored  for  a  problem  like  the  present 
one. £21“221  The  first  part  of  that  problem  would  be  to  determine  a  steady-state 
flow  pattern  for  a  fluid  with  a  distributed  heat  source  within  a  right 
circular  cylinder  with  Its  axis  aligned  at  some  angle  to  the  vertical.  For 
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the  case  of  a  horizontal  cylinder  of  Infinitely  small  radius,  the  flow  pattern 
has  been  calculated  by  Y1 h.  C23-243 

Unfortunately,  however,  that  solution  is  not  of  great  value  for  the 
present  problem  because  the  size  of  the  cylinder  radius  Is  a  critical  parameter. 
Nevertheless,  It  Is  expected  that  Ylh's  solution  will  assist  In  obtaining  the 
asymptotic  form  of  the  steady-state  flow  at  large  distances  from  the  laser 
beam  cylinder,  fnce  that  time- Independent  flow  pattern  has  been  determined, 
the  linearization  of  the  hydrodynamic  equations  for  perturbations  from  the 
flow  pattern  will  lead  to  an  eigenvalue  problem,  which  eventually  will  yield 
a  critical  value  for  W.  Ostrach^25^  suggests  that  the  eigenvalue  problem  can 
be  bypassed  as  the  stability  of  fully  developed  natural  convection  flows  can 
be  found  by  using  the  appropriate  velocity  profile  In  the  classical  theory  of 
hydrodynamic  stability.  This  assertion  rests  upon  his  analysis  of  the 
stability  of  free  convection  above  a  flat  heated  plate,  where  Instability 
first  appears  for  a  Reynolds  number  of  283. 

Above  the  threshold  for  beam  Induced  turbulence,  governed  by  Wcr, 
general  arguments^26^  lead  to  a  size  for  the  smallest  eddies,  a(W  /W)3^. 

V  ■ 

3 

For  a  1-meter  beam.  If  Wcr  should  be  about  10  ,  then  the  eddies  might  have 

2 

sizes  as  small  as  7  cm  for  a  power  flux  of  1  kW/cm  .  The  associated  density 
fluctuation  would  then  be  expected  to  result  In  considerably  Increased 
scattering  of  the  beam. 

The  arguments  presented  here  show  that  there  are  substantially  more 
Important  sources  of  Instability  In  the  laser-fluid  system  than  those  discussed 
In  earlier  linearized  analyses.  It  Is  felt  that  these  fluid  Instabilities  will 
be  enhanced  by  their  Interaction  with  the  scattering  of  the  laser  beam,  because 
of  the  general  result  that  Instabilities  In  fluids  result  If  the  heating  of 
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the  fluid  is  greater  in  those  regions  where  the  density  of  the  fluid  is 
greater. 

At  the  present  time  we  can  only  outline  the  general  nature  of  the 
effects  to  be  exnected  above  a  critical  power  level.  Much  additional  work 
clearly  needs  to  ije  done,  both  of  an  experimental  and  theoretical  nature. 
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PART  VI 

LASER  RADIATION  IN  A  MOVING  MEDIUM 


In  order  to  consider  the  effects  of  forced  fluid  flow  on  a  laser  beam, 
the  equations  for  the  laser-fluid  system  will  be  analyzed  for  a  two-dimensional 
steady  state  situation.  For  simplicity,  consider  a  uniform  slab-shaped  laser 
beam  of  width,a, in  the  y-direction,  of  infinite  extent  in  the  x-direction, 
and  propagating  in  the  z-direction.  If  there  is  a  transverse  motion  of  the 
medium  in  the  y-direction,  due  to  wind  or  to  convection  above  a  heated  surface, 
an  asymmetric  distribution  of  temperature  and  density  will  be  generated  in 
the  beam  region.  This  asymmetry  leads  to  a  curving  of  the  laser  beam  toward 
the  source  of  the  fluid  motion  as  shown  in  Figure  [18]. 

This  self-curving  or  "wind  prism"  effect  can  be  explained  by  the  follow¬ 
ing  simple  physical  argument.  Consider  the  case  of  fast  flow,  where  flow 
transit  time  through  the  beam  cross  section  is  much  less  than  thermal  conduction 
times.  [The  flow  is  not  presumed,  however,  to  move  at  speeds  comparable  to  the 
sound  velocity  in  the  medium.]  As  the  fluid  passes  the  beam  cross  section,  it 
is  heated  and  becomes  less  dense,  so  that  the  density  at  y  ■  a  is  less  than  at 
y  *  0.  Thus,  the  index  of  refraction  at  y  =  a  is  less  than  at  y  ■  0,  so  that 
the  beam  is  effectively  moving  through  a  prism  and  will  bend  toward  the  negative 
y-direction.  As  this  wind  velocity  decreases  to  zero,  the  wind  prism  goes 
away,  to  be  replaced  by  an  effective  convex  thermal  lens.  In  the  absence  of 
wind,  thermal  conduction  will  lead  to  a  temperature  distribution  symmetric 
about  y  =  j.  The  hottest  place  will  be  in  the  center  of  the  slab  and  the 
density  distribution  will  have  a  minimum  there.  Thus,  the  index  of  refraction 
has  a  minimum  at  y  *  j  and  the  beam  is  effectively  moving  through  a  defocusing 
lens.  The  beam  will,  therefore,  tend  to  spread  out  or  "bloom"  as  indicated  in 
Figure  [19]. 
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Akhmanov,  et  al.,^8^  and  other  groups^6'17'28”31  ^  have  predicted 
and  observed  the  phenomena  of  thermal  blooming  and  self-curving  In  the  presence 
of  a  wind.  These  two  effects  will  be  considered  simultaneously  here  for  the 
simple  slab  geometry  and  a  formula  for  the  development  of  the  beam  profile 
will  be  given  which  included  both  types  of  behavior.  The  slab  geometry  will 
overemphasize  the  self-curving  effect  slightly,  compared  to  a  cylindrical 
beam,  but  qualitative  features  and  order  of  magnitude  will  be  correct.  The 
results  of  the  present  model  will  be  compared  with  a  cylindrical  beam  experiment 
at  the  end  of  this  section  and  reasonable  agreement  will  be  demonstrated.  In 
fact,  this  comparison  will  serve  to  illustrate  the  importance  of  thermal 
conduction  in  the  description  of  a  "blooming  wind  prism." 


In  order  to  relate  the  two  effects  and  to  make  quantitative  estimates 
of  the  curvature  of  the  beam,  equations  (25a)  -  (29)  for  the  laser-fluid 
system  will  be  specialized  to  this  case.  A  cartesian  geometry  is  used  with 
the  fluid  variables  T  and  p  being  functions  of  y  only,  and  the  fluid  velocity 
v  =  v(y)ey  being  a  function  of  y  only  and  directed  in  the  y-direction.  These 
variables  are  also  taken  to  be  time  independent  since  only  a  steady-state 
analysis  Is  to  be  attempted.  Similarly,  the  electric  field  will  be  taken  to 
be  of  the  form 


E  =  E(y)e  e 


i(u>Lt-kLz) 


01  7 

1Z 


+  c.c. 


(223) 


where  e  is  a  polarization  unit  vector,  and  k^  are  the  frequency  and  wave 
number  of  the  laser  radiation,  and  the  factor 


.  2  z 


included  to  allow  for  damping  consistent  with  the  damping  term  in  the  wave 
equation,  (25a).  Substituting  (223)  into  (25a),  an  equation  for  E(y)  is 
obtained: 
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where,  by  (25b), 


A 

37 


(k* '  7 c) 


(224) 


e  *  +  &2  ^ 


(225) 


and  Is  easily  computed  from  (223): 


7?  -  1  rZi 


E‘  =  ?  E4(y)  e-az  (226) 

To  further  describe  e  as  a  function  of  p  and  T,  the  Clausius-Mossotti  relation 
will  be  used  as  shown  in  (25c): 


/3e\  ~  (e-1)  (e+2)  /9e\  ^ 

Mt  *  ‘  ~ 


(227) 


Since,  for  air  and  many  other  fluids,  e  is  independent  of  T  for  reasonable 
temperatures,  (  227)  can  be  integrated  at  once  to  give 


e  = 


(eo*2>P0  * 


(228) 


where  en  =  [e.  +  e,E^]  and  pn  are  the  values  of  e  and  p  evaluated  at 
0  L  £  y=0  0 

y  *  0.  This  relation,  also,  will  be  called  the  Clausius-Mossotti  relation. 


Now,  if  p  can  be  found  as  a  function  of  y,  then  (228)  and  (224)  can  be 
used  to  compute  E(y)  and  the  ray  methods  of  geometrical  optics  can  be  used  to 
describe  the  deflection  of  the  beam.  From  the  wave  equation,  one  obtains  the 
following  relation  between  the  gradient  of  e  and  the  radius  of  curvature  Rc  of 
the  laser  beam 


vi"  ■  7i  * 


(229) 
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where  V^n  Is  the  transverse  component  of  the  gradient  of  the  index  of 
refraction.  The  deflection  of  the  laser  beam,  therefore,  can  be  described  by 
the  relation 


(230) 


where  ip  is  the  scattering  or  deflection  angle  shown  in  Figure  [18]. 


In  order,  then,  to  compute  p,  the  Navier  -Stokes  and  heat  transfer 
(conservation  of  energy)  equations  (26a)  -  (27b)  must  be  solved  along  with  the 
fluid  continuity  equation  (28)  and  the  equation  of  state  of  the  fluid  (29). 
Specializing  these  equations  to  the  present  situation,  one  gets: 


NAVIER- STOKES 


u  dv 
pv  37 


dP  .  .  ,,  d2v  .  (e-l)(e+2)  ll  dE2  .  e-1  7  dp 

'  37  *  (2n+n  1  ^  3^ - ‘[TF'TIdf 


1 


(23l) 


HEAT  TRANSFER 


v 


dT  £0 
dy  '  p8 


v 


dp  _  2n+n‘ 
dy  “  pCv 


<  d2T 
pCv  dy2 


+ 


ac/e  i-2 

pC  L 
H  v 


+ 


FLUID  CONTINUITY 

pv  =  pQvQ  =  constant 


(232) 


(233) 


EQUATION  OF  STATE,  IDEAL  GAS 

P  =  RpT  (234) 

where  pQ  and  vQ  are  the  density  and  wind  velocity  .at  y  -  0  and  R  is  the  gas 
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constant  In  units  of  ergs  per  gram  per  degree.  To  obtain  (231),  the  Clatilus- 
Mossottl  relation  was  used  and  the  dependence  on  gravity  was  dropped.  For 
transverse  winds,  the  force  of  gravity  would  act  In  the  x-directlon  and  would 
destroy  the  simplicity  of  the  present  model.  For  "winds"  caused  by  rising 
columns  of  air  over  heated  surfaces,  however,  gravity  could  act  in  the  y- 
dlrectlon  and  would  cause  no  essential  complication  in  the  present  model.  For 
simplicity,  the  gravity  term  will  simply  be  dropped  in  the  analysis  given  here. 

Before  proceeding,  several  additional  simplifications  will  be  made  In 
(231)  and  (232).  The  terms  Involving  the  viscosity  coefficients,  (2n+n')» 
are  very  small  and  will  be  dropped.  Furthermore,  the  terms  involving  deriva¬ 
tives  of  the  thermal  conductivity  k  are  small  and  will  also  be  dropped.  The 
continuity  equation,  ( 233) •  and  the  gas  law,  (234).  can  be  used  to  eliminate 
P  and  v  from  (231)  and  (232).  Thus,  the  density  and  temperature  are  described 
by  the  two  equations: 


(p0v»> 


o1  dp  .  R  d(pT)  (e-l)(e+2)  fl  d?  e-1  7  dpi  ,  , 

dy  ciy  T  \J  dy^  ^  P  <tyj 


and 


vo3y  ‘  v 


T  dp  k  d^T  ,  qc/c  t: 

o  p  dy  p  C  .  2  p  C 

M  3  Ho  v  dy  Ho  v 


(236) 


where,  in  order  to  obtain  (236),  the  relation 


(237) 


following  from  (234),  was  used. 

Inspecting  (236)  for  the  characteristic  times  for  thermal  changes. 
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one  observes  two  characteristic  times:  tc>  the  thermal  conductivity  time, 
and  Tp,  the  flow  transit  time  through  the  beam  cross  section.  These  times 
are  conveniently  taken  to  be: 


C  a2p 
_  v  o 

C  " 


(238) 


and 


(239) 


For  convenience,  these  time  scales  will  be  compared  to  the  much  shorter  sound 
transit  time: 


T 


S 


(240) 


where  v$  is  the  isentropic  velocity  of  sound  at  temperature  TQ.  It  will  also 
be  convenient  to  Introduce  the  following  dimensionless  parameters: 

w  =  * 
w  -  a 

0  =  temperature  relative  to  temperature  at 

'o  y  =  0 

g  =  =  density  relative  to  density  at  y  =  0 

po 

P  -p 

f  =  1-g  =  —  =  relative  density  decrement 

po 

Using  this  notation,  (235)  and  (236)  can  be  written  in  the  form: 

1  dE2  .  €-1  E2  dol 

iair  +  -rrd;J  (245) 

and 


1 

2 

Tp 


dJi> 

dw 


1 


YT. 


d(qe)  _  (e-1 ) (e+2) 
2  dw 


3o_a 


(241) 

(242) 

(243) 

(244) 
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(246) 


]_  d0  (v-1  )  0  d£  _  1__  aQ  + 

tf  a?  -  tf  g  a? '  tc  ^7 


ac/e 

c  d  T 
vMo  o 


? 


Now,  the  general  strategy  will  be  to  solve  (245),  (246)  and  (224) 

by  Iteration,  but  only  one  Iteration  will  be  made.  A  simple  shape  will  be 

presumed  for  E(y)  and  then  (245)  and  (246)  will  be  solved  for  0  and  g.  Then, 

In  principle,  one  could  solve  (224)  for  a  new  E(y)  and  continue.  However, 

the  first  Iteration  for  0  and  g  allows  evaluation  of  Ve,  so  that  (230)  can 

be  integrated  to  produce  an  estimate  of  the  beam  deflection.  The  initial 
f 

shape  presumed  for  E(y)  is 


E(y)  = 


=  constant,  for  0  <_  y  <_  a 


(247) 


(0  otherwise 

"7 

Then,  Inside  the  region  0  <_y  <  a;  i.e. ,  0  <  w  <  1,  the  derivative  of  E 
with  respect  to  w  is  zero  in  (245).  In  both  (245)  and  (246),  the  replacement 


7?  _  1  f2  -az 
E  -  ?  Eo  e 


\  _A e 


•az 


7e~"  c 
o 


(248) 


can  be  made,  wnere  IL  is  the  power  per  unit  area  in  the  incident  laser  beam. 
For  convenience,  now,  the  following  additional  time  scale  is  introduced. 


tA  = 


C  p  T  C  p  T 

V *0  O  _  v  0  0 


C  p  T 
vHo  o 


acv^  E^ 


al, 


al^e 


-az 


(249) 


The  time.T^,  Is  characteristic  of  the  rate  of  absorption  of  heat  in  the  fluid. 


+ 

Such  a  flat  profile  will  minimize  the  thermal  blooming  because  none  of  the 
rays  will  cross  in  the  first  iteration.  Subsequent  iterations  would  produce 
crossing  rays  which  would  keep  the  iteration  procedure  from  converging. 
Crossing  rays  are,  of  course,  not  allowed  by  the  wave  equation. 
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Now,  (245)  can  be  integrated  in  the  following  way 


-7(98-1)  =  - 7  I"  (e'U(ct2)  (e-l)  g*  dw  = 

tf  9  n2  3-^%a2Jo  39 


[wd£(e.i)dadws 

3/T  cp  a2  Jo  ^  * 


0  ^0 


L= _  fW  dkUl  da  dw  . 

“  Cn  a2  Jn  d9  dw 


6*^o  Cpoa  0 


6t/e~  cp  a 
0  0 


From  (228),  however,  one  gets 


L - jCte-l)2]’ 


(e-l) 


-  (  38q3  V- 
'  iFw 


2  .  I  "JUT.  [3B0g(l+80g)]2  -  9B2  (g2  +  2B0g3) 


where 


e-l 

B0  =  —jg  ~  very  small  for  air 


Combining  (250)  and  (251)  and  solving  for  0  produces 

where  is  still  another  characteristic  time: 

2  _  cpoa2 


D"  K\ 


(250) 


(251) 


(252) 


(253) 


(254) 


This  quantity,  Tp,  characterizes  the  direct  interaction  of  the  electromagnetic 
field  with  the  fluid  and  arises  from  the  electromagnetic  stress  term  in  the 
Navier-Stokes  equation. 
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Equation  (253)  gives  the  temperature  as  a  function  of  the  density  and 
will  be  useful  In  Integrating  (246).  If  (253)  Is  used  to  express  the  ratio 

g 

-as  a  function  of  g,  then  a  first  integral  of  (246)  can  be  obtained  immediately. 
Then  (253)  can  be  used  to  eliminate  9  in  this  first  integral  and  a  first  order 
differential  equation  for  g  results.  Continuation  along  this  line  is  largely 
superfluous,  however,  and  leads  to  implicit  transcendental  expressions 
connecting  p  and  w.  It  is  useful, and  far  simpler,  to  agree  to  study  the 
situation  for  the  case  in  which  the  density  change  is  small  from  y  *  0  to 
y  ■  a. 


For  the  case  that  g  2  1,  the  relative  density  decrement,  f,  defined  in 
(244),  is  very  small  and  it  is  useful  to  write  equations  (253)  and  (246)  in 
terms  of  f  and  work  to  first  order  in  the  density  decrement.  Thus  (253) 
becomes 

0  ~  1  +  Qf  (255) 

where  2  2 

q  =  1-Y  (£)  .  2y  Q  <H3B0)  .  (255) 

Then,  using 

i:  =  !  +  (Q+,)  f>  (257) 

equation  (246)  can  be  integrated  immediately: 

Ldl.iiiiilif  s.H_+A]  (258) 

where  A1  is  a  constant  of  integration.  Integrating  again  and  taking  the 
boundary  conditions  f(0)  =  0  and  p(l)  =  p^,  where  p-|  is  the  density  of  the 
fluid  at  y  =  a,  produces 
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f  = 


Vfl  tf 

c0  YTA(l-«) 


exw-  1 


TpW 


,x  .  1  +  YTaO-o) 


(259) 


where 


and 


X  =  Y  “  1-° 


(1+3B0) 


tf  1-y a 


(260) 


(261) 


Then,  combining  (251)  and  (259)  the  square  root  of  the  permittivity  is  approx¬ 
imately  given  by: 

il+B_ 


Jz  ~  1  + 


I  Bo  (,+2Bo) 


1+2B. 


(262  ) 


and 


n  r.  &Jz  df  dw  _ 

VE  ‘"df  did7' 


Since  this  derivation  has  proceeded  under  the  assumption  that  the  density 

Vpl 


decrement,  f,  is  small,  the  first  order  term  involving 


is  small  compared 


to  1  and  will  be  dropped  as  will  terms  of  order  BQ  when  compared  to  1.  It  is 
convenient  now  to  introduce  the  notation: 


tF  tF 


— -  =  p  1L.  =  fem 
YtA  YTC 


where  the  dimensionless  power  fi  is  given  by 


and 


(264) 


fi  =  fie 


az 


(  265) 


(263) 
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The  deflection  of  the  laser  beam  can  now  be  given  in  the  form 


dili  •  (eo_1)^ 

Hz  ~  S  TT^y 


1  e 


X*  1 


x  ex-l 


(266) 


In  order  to  analyze  this  expression  for  the  curvature  of  the  laser  beam,  it  is 
convenient  to  rewrite  o  In  the  form 

‘2.  -i 


.(if -A 

\vs  / 


(267) 


0  C  V 

0  Ho  s 


For  propagation  through  air  at  10°C,  the  two  terms  on  the  right  hand 

-3  2 

side  of  (267)  are  both  of  order  10  when  vQ  ~5  (km/hr)  and  IL  *  1  (MW/cm  ), 

so  o  and  yo  can  be  considered  small  compared  to  1  for  velocities  less  than 

50  (km/hr)  and  powers  less  than  100  (MW/cm  ).  From  equations  (261),  (238), 

and  (239),  one  has 


X  * 


/  l-o\  .  /  l-o\  *V>o 
h^j  Tf  "  y’Ro/  — IT  "o 


(268) 


so  that  for  wide  beams  or  fast  winds,  the  conduction  time  is  very  long,  com¬ 
pared  to  the  flow  time  and  X  will  be  very  large  compared  to  1.  Under  these 
conditions,  the  ratio  of  exponentials  in  (266)  tends  to  be  smaller  than  j  , 
so  that  most  of  the  rays  bend  into  the  wind.  For  large  enough  y,  however,  the 
ratio  of  exponentials  is  greater  than  ,  so  that  some  of  the  rays  have 
negative  curvature.  Thus  the  downwind  edge  of  the  beam  bends  in  the  direction 
of  the  wind.  Therefore,  pure  beam  curving  with  no  thermal  blooming  will  never 
occur. 


Physically,  this  is  easy  to  understand.  The  hottest  place  and,  hence, 
the  minimum  of  the  index  of  refraction,  must  lie  inside  the  beam  region.  On 
the  two  sides  of  this  minimum,  the  rays  will  bend  in  opposite  directions.  The 
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faster  the  wind,  the  larger  X  will  be  and,  hence,  the  larger  y  must  be  before 
the  curvature  changes  sign.  In  other  words,  the  wind  can  blow  the  hottest 
spot  away  from  the  center  of  the  beam,  but  it  can  not  blow  the  hottest  spot 
entirely  out  of  the  beam.  The  position  of  the  minimum  of  the  index  of 
refraction  is  easily  obtained  by  determining  the  value  of  y  for  which  the 


right  hand  side  of  (266)  vanishes: 

yo  _  ,  log  X 
a  X 

i  mg  llz£h 

X 

(269a) 

«1  .  iQf  A 

A 

for  large  X 

(269b) 

_,1  .  X 

“2  +  2T 

for  small  X 

(269c) 

In  particular,  for  very  small  X  the  minimum  is  at  the  center  of  the  beam, 
whereas,  for  large  X,  the  minimum  is  very  near  the  downwind  edge  of  the  beam. 


As  the  wind  velocity  goes  to  zero,  equation  (266)  reduces  to 

_  (yO?  r, 


dA  .  _____ 
dzjv  =0  2a  0-Yo) 


M 


(270) 


which  is  symmetrical  about  the  center  of  the  beam  as  would  be  expected  for  pure 
thermal  blooming.  The  rays  at  the  edge  of  the  beam  have  the  maximum  curvature 


max 

d^ 

_  1 

y 

dz 

|CM 

II 

O 

II 
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> 

2a  (1-yo) 

(271 ) 


In  terms  of  this  maximum  curvature  for  pure  thermal  blooming,  (266)  can  be 
written  in  the  form 


1  ® 

x  " 


ex-l 


(272) 


Clearly  £  is  a  monotonic  function  of  y,  giving  the  largest  positive  curvature 
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at  y  =  0: 


5 


max 


«(1  -  -g-)  for  small  X 


h 


for  large  X 


(273a) 

(273b) 

(273c) 


As  may  be  seen  from  (273a) ,  £max  decreases  monotonlcally  with  the  wind  velocity 
and  Is  less  than  5B  If  there  Is  any  wind  at  all.  Thus*  the  wind  blows  the 
blooming  rays  back  Into  the  beam  on  the  upwind  side  of  the  beam. 


The  minimum  curvature;  that  Is,  the  largest  negative  curvature,  occurs 


—  <’  *y>  h 
0  -  T>  h 


h 

(274a) 

for  small  X 

(274b) 

for  large  X 

(274c) 

Thus,  It  is  clear  that  Increases  monotonlcally  with  the  wind  velocity 

and  Is  greater  than  5g  if  there  Is  any  wind  at  all.  The  effect  of  the  wind  on 
the  downwind  side  of  the  beam  Is  to  blow  the  blooming  rays  farther  out  of  the 
beam. 


Comparing  (273)  and  (274),  one  finds  that 


W  -  cmi„  ■  W  ♦ 


*m1n 


=  25, 


(275) 


regardless  of  the  wind  velocity.  Thus  the  wind  has  no  effect  on  the  size  of 
the  spread  of  curvatures  created  by  the  thermal  blooming.  For  v  =  0,  this 
spread  Is  from  5g  to  -5B,  whereas  for  extremely  fast  winds,  the  spread  Is 
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essentially  from  0  to  -2£g.  Thus,  the  curvature  spread  can  not  be  decreased 
by  a  wind,  but  it  can  be  translated  by,  at  most,  half  nf  its  extent.  There 
is,  however,  a  very  dramatic  shift  in  the  fractions  of  the  power  being 
scattered  into  the  various  curvatures.  This  arises  because,  for  large  X, 
the  distribution  shown  in  (272)  is  very  flat  as  a  function  of  0  for  y  <  y  . 
For  X  =  100,  for  example, 


for 


0.96. 


Then,  for  y  >  yQ,  the  curvature  becomes  negative  and  falls  rapidly  to 


-  -  0.99 


at  y  =  a.  The  net  effect,  therefore,  Is  that  about  96%  of  the  power  goes 
off  with  curvature  £  ~  Cg/50  and  the  remaining  4%  blooms  strongly  out  of  the 
beam.  For  a  non-uniform  beam  profile,  the  conclusions  are  qualitatively  the 
same,  but  the  distribution  is  no  longer  monotonic  as  in  (272),  so  that  some  of 
the  rays  will  cross. 


One  of  many  ways  to  exhibit  the  nature  of  the  scattering  distribution 

i 

given  in  (272)  Is  to  show  the  scattering  cross  section  for  scattering  of  power 
into  the  various  curvatures.  Thus,  one  can  define  the  cross  section 


o(£)d5  =  ap(c)dc 


laser  energy  flux  scattered  into  range  dg 
laser  energy  flux  density 


ad£ 


(276) 


The  quantity  p(5)  is  the  fraction  of  the  laser  power  scattered  into  curvature 


S.  The  distribution  p{£)  is  normalized  to  unity: 
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’max 


P(5)d£  = 


E  . 
smin 


’max  di  =  ^  ^  2Zb  -  XSmin  _  1 

X 


2Kb-XK 


2h  -  X?max 


j  log  e  =  1 


( 2/7) 


’min 


where  £L.  and  are  given  in  (273a)  and  (274a).  The  fraction  of  the  power 
min  max  a  r 

being  scattered  into  negative  curvatures  is 


(278) 


The  distribution  p(£)  is  sketched  In  Figure  [20]  for  small  wind  speeds. 
The  extreme  values  of  the  monotonic  distribution  are  given  by 
by 


ex-l 


h  PUmax}  =  ~2T 


and  (  279) 

P^miV  =  1  2X6  * 

•r 

As  shown  in  the  figure,  the  distribution  is  sharply  enhanced  at  ?  ,  even 

for  X  =  3.  Because  of  this  sharp  freaking  of  pU)  for  even  very  slow  wind 
speeds,  the  thermal  lens-wind  prism  effect  may  be  pictured  as  in  Figure  [21]. 
At  each  point  In  the  beam  path  the  beam  Is  essentially  passing  through  a  wind 
prism  with  a  small  strongly  defocusing  thermal  lens  perched  on  the  downwind 
vertex  of  the  prism.  As  X  increases,  the  half-angle  of  the  prism  decreases 
rapidly  and  the  blooming  lens  Increases  its  defocusing  power,  but  shrinks 
rapidly  in  size  so  that  only  an  Infinitesimal  portion  of  the  downwind  edge  of 
the  beam  is  affected.  As  the  wind  velocity  Increases,  both  the  curving  and 
the  blooming  effects  decrease  as  indicated  In  equations  (273c)  and  (278). 
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one  has,  from  (268): 


1400  a 


2  vo 


"1 


*  1400 


so  that  at  a  given  point  in  the  beam,  about  0.5%  of  the  beam  is  affected  by 
the  thermal  lens.  The  dimensionless  power  shown  in  (265)  is 


P-  0. 


043  a£  *  0.043 


-6  -1  ^ 

where  the  values:  a*  3  x  10  cm  ,  k  *  2.5  x  10  ergs/ (deg  sec  cm),  and 
Tq  *  283°K  =  10°C  have  been  used.  The  blooming  edge  curvature  is  then  given 
by  (271): 

e  ”1  -4 

?b  ~  (°-°43  a2>  ~  ~~H  1‘° . {°*043  a)  =  6‘04  x  10"6  a  ~ 


s  6.04  x  10 
so  that,  from  (273c), 


-6 


Cmax  =  (8-64  x  10'9)  a  -JE:  8.64  x  10’9  . 


10 

7 


km 


Thus,  after  the  beam  has  propagated  1  Ion,  it  will  have  an  inclination  of 

-4 

the  order  of  9  x  10  radians  ~  0.05  deg.  For  comparison,  without  the  wind 
the  pure  thermal  blooming  would  have  produced  a  70°  divergence  of  the  beam! 
The  ratio  of  \|/„w  to  is 

UkJX  D 

8.- 8.4,  x JO  9  .  1>43  x  10-3 
6.04  x  10'° 

so  that  the  wind  can  play  a  significant  role  in  holding  the  beam  together. 


The  fraction  of  the  power  lost  from  this  beam  in  1  km  is  roughly 

i ... 

T. 


V™  -  /,  -OZx  ,  i  , 

—  Cl  -  e  )  +  j  log  — — 


(1  -  e'0,3)  +  ^  log  [6.04  x  104]  =  0.258  +  0 


.008 
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*  “CtZ 

so  that  the  major  loss  of  power  arises  from  the  damping  factor  e  in  the 
"prtSfance  of  the  10  j^-  wind. 

It  is  possible  to  compare  the  predictions  of  this  section  with 

experiment,  using  the  nice  results  of  Smith  and  Gebhardt.^J  The  beam  of 

a  20  watt  cw  C02  laser  operating  in  the  fundamental  TEMQ0  mode  was  passed 

through  an  enclosed,  recirculating  wind  tunnel  with  a  100  cm  optical  path 

length  and  a  variable  wind  speed  from  0  to  500  cm/sec.  The  absorption  of 

the  air  was  increased  more  than  a  thousandfold  and  was  controlled  by  the 

addition  of  varying  amounts  of  freon-12  (CC1 2^2^  *  ^  vert*cal  opt^31  path, 

transverse  to  the  wind  direction,  was  used  to  minimize  the  gravitational 

convection  effects.  The  -  beam  diameter  was  0.4  cm  at  the  tunnel  entrance 

e 

window  and,  in  the  absence  of  thermal  distortion,  the  ^-diameter  was  0.6 
at  the  scanning  detector  which  was  located  20  cm  from  the  wind  tunnel  axis 
window.  Of  course,  the  experiment  was  performed  with  a  cylindrical  beam, 
but  one  might  expect  to  get  semiquantitative  agreement  from  the  present  slab 
beam  model,  particularly  if  one  concentrates  on  the  shift  of  the  beam  peak 
into  the  wind. 


To  consider  the  predictions  of  the  slab  beam  model  for  the  experiment 
of  Smith  and  Gebhardt,  one  can  use  equation  (272).  A  ray  initially  at  height 
y  will  have  the  trajectory  Y(y,z),  where 


37  =  -tan  ij )  ~  -ip. 

Differentiating  again  and  using  (272).  one  gets 

d^Y  d£  5  i  -f  (y)  e“az  for  z  <  L 
dz  '  0  for  z  >  L 


(280) 


(281) 
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where  L  is  the  tunnel  length  (100  cm)  and 


(282) 


Integrating  (281),  one  has 


Hi¬ 


tt 


dY  _  f(y)  ,  J(1  -  e_aZ)  for  2  <  L 


(1  -  e’aL)  for  z  >  L 


(283) 


Another  integration  produces: 

Y(y,z)  =  y  -  [*(1  +  aL) 

cl  L 


-aL 


■i 


-  1  +  (1  -  e"aL)  azl  for  z  >  L  (284) 


In  particular,  for  z  =  ^  L  =  100  cm  +  20  cm  for  L  =  100  cm,  the  height  of  the 
ray  at  the  detector  will  be 


YDet  _  Y(y’  5  L)  _  y  2*’BaL  f  n  [l  e*9  1 

—  =  — 1-  -a--a-2T9(aL)b‘7TJ 


(285) 


where 


9(ol)  *  [l  ♦  (i  •  i)  (l  •  *'“L)]  • 


The  corresponding  intensity  at  z  «  L  is 
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(286) 


The  peak  in  this  intensity  distribution  occurs  for  y  -  0  and  its  location  is 
determined  from  (285)  to  be 


L>. 

a 


2^b°(L 

a2a 


g(aL) 


1_  J-l 
X  eX-lJ 


(287) 
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This  relation  is  conveniently  put  in  the  form 


G(al) 


Ie 

a 


where 


G(aL)  = 


(288) 


(289) 


In  order  to  compare  the  shift  of  the  peak  with  the  results  of  Smith  and 

Gebhardt  It  will  be  necessary  to  put  numbers  into  (288).  The  numerical  values 

for  the  various  parameters,  appropriate  for  TQ  *  10°C  are  shown  in  Part  II  of 

this  report.  The  width,  a,  of  the  slab  beam  will  be  taken  to  be  0.4  cm,  the 

^  diameter  of  the  experimental  beam.  The  20  watt  beam  used  would  have  a 

"uniform"  transverse  linear  density  of  20  watt/0.4  cm  so  that  IL  will  be  taken 
o 

to  be  50  watts/cm  .  Then,  in  cgs  units,  one  has 


Thus,  for  the  experiment  under  consideration,  (288)  can  be  written 


-G(otL)  =  (9.95) 

a 


^65 


»$n>  v. 


- 1  j 


(290) 


In  Figure  [22]  the  function  G(aL)  is  plotted  versus  (aL)  and,  also, 
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-  6(ocL )  ,  the  RHS  of  (290)  is  plotted  versus  [■—-)•  Data  given  by  Smith  and 

o 

Gebhardt  for  aL  =  0.14,  0.53,  and  1.0  are  also  indicated  on  the  figure.  One 
notices  that  the  data  for  various  aL  does  not  all  lie  on  the  same  "universal" 
curve  as  suggested  by  (290).  Thus,  as  expected,  there  is  some  difference  in 
a  slab  beam  and  a  cylindrical  beam.  The  important  feature  to  be  noticed, 
however,  is  that  the  order  of  magnitude  of  the  predicted  effect  is  correct 
and  the  "saturation"  of  the  prism  effect  encountered  by  Smith  and  Gebhardt  is 
predicted  by  (290).  It  would  certainly  appear  that  thermal  conduction  effects 
can  account  for  the  saturation  of  the  deflection  and  that  the  other  explanations 
suggested  by  Smith  and  Gebhardt  are,  perhaps,  unnecessary. 

In  Figure  [23]  the  intensity  distribution  of  the  laser  pulse  at  the 
detector,  determined  by  equation  (286),  is  shown  for  vQ  =  2  ,  10  ,  and 

for  vQ  extremely  large.  One  notices  here  the  rapid  transition  from  the 
uniform  spreading,  due  to  thermal  blooming,  to  much  narrower  distributions, 
peaked  on  the  upwind  side,  as  the  wind  velocity  increases.  The  qualitative 
features  of  these  curves  were,  of  course,  already  illustrated  in  Figure  [20]. 

This  analysis  has  been  undertaken  in  order  to  estimate  the  relative 
importance  of  thermal  blooming  and  wind  curving.  Certain  qualitative  features 
have  been  discussed  and  perhaps  some  insight  has  been  achieved.  The  approx¬ 
imations  made  have  been  rather  strong  and  cylindrical  beams  and  non-uniform 
beam  profiles  have  not  been  properly  considered  here.  These  matters  are 
currently  under  study  and  results  will  be  presented  at  a  later  time.  The 
approach  to  steady  state  and  the  usefulness  of  working  to  first  order  in  the 
density  changes  is  also  under  study.  There  are  difficulties  involved  in 
asking  for  steady  state  solutions ^4nd  the  effect  of  a  wind  on  such  matters 
must  be  analyzed  carefully.  It  does,  however,  appear  that  blooming  in  the 
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absence  of  a  wind  is  a  quite  serious  problem  if  one  wishes  the  beam  to 
retain  its  character  over  long  distances.  A  fast  wind  can,  apparently,  help 
considerably  to  alleviate  the  thermal  blooming. 
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APPENDIX  A 

Computer  Program  for  Power  Threshold  for  Instabilities 

One  can  show  from  equation  (68)  that  the  threshold  power  for 
an  instability  li*»s  near  k=0.01cm-^.  A  search  was  conducted  for 
various  v  with  k  fixed  at  O.Olcm"^.  The  power  threshold  was  then 
located  at  y=(8.45xl0“2)  *k  and  found  to  be  Ptj1=312ergs/sec. 

In  this  computation  the  following  notations  have  been  intro¬ 
duced  into  the  computer  program: 

S  :  k 

R  :  k3 

B (6)  :  Coefficient  of  the  5th  order  polynomial  obtained  from  (68). 

ROOT (5)  :  The  roots  of  the  dispersion  relation. 

DET  :  Decrement  of  the  power. 

DCROOT  :  Subprogram  for  the  calculation  of  the  roots  of  a 
polynomial. 

Following  is  the  computer  program  for  this  computation. 
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wcimm: ppiir  options  imaimi: 

net  tw.nw.i  .I'cmfi  ,ri  ,F?,r^,r',,F0,P6,p?,r*j,r‘>«  cplxiiai  ; 

M»i.u  vs*sortu.i7si'm: 

nr. ii oim  iMinc  (a.u.aci. ,Rom,  ir.iiK.i'Kr.i  Ni  i 

/♦THIS  1 A  A  SURROIJT  1  Ml-  FUR 

THi:  ROOTS  OF  A  COMPLEX  MOLYNOMUl  OF  OKCRFF  LESS  TUAN  NXY  WITH  16 
OECIMAL  OIGITS  ACCURACY  IN  STORAGE*/ 

MXf-61 

ML  I  A I  6  (  ,70,21 ,72, 7x,H?,H3tXL?tXL3,n?lG?,F7Q(F7L,F72,Kmim< 

HUF  f  41 ,  RE  OfNI  61  I  Cl'ix  1(6(1 

IF  All  I  *  .0  TUFIJ  no:  RIM  rat  ISYSPRINTI  list-  (•IFA01IIG  COEFFICIENT'  is  IP 
RO  |M  nCPUI’T  *  I ;  FXIT:  FNO;  IF  N  >*  (IXY  I  14  <  t  THIN  00*.  PUT  r  IlK  ISYSPRINTI 

list  ( •  pucp ec  nr  polynomial  in  ocrout  is  mot  retwcen  o  and  nxym; 

rx(T;  END; 

NR»ns  K* Mf  l :  RI'rGCN  *  Al 

Si:  ZO  — i.s  7  l»l . ;  7  ?*.n;  XI.?*-. 5}  02*. S«  00  l*M  TO  l  RY  -1: 

IF  All  (  * .  0  TIIF.N  on:  MR  *  NR  ♦  I  ;  ROOT  (UR  (=  .0;  END;  FI  \E  GO  TO  S7  ;  ENni 
ST:NN*l:  (F  101  2  THEN  GO  TO  S?5:  FI  SC  IF  NN*2  THIN  DO:  7  3=- AI2  |  /  Al  \  I  J 

NR*NRf  1  *,  P00TIMN*;V,  GO  TO  S2D5  I  MO  ■  F  70*A»  NNI -A  IUU-1  I  ♦  Al  NN-  ?»  5 
F 1 1  ■  A I WM I  ♦  M IIN-  1 I ♦ A (NN -2 I  I  F2?  =  A(liNI;  IIN*NN-1  i 
S3!  G2*r77*l  XI  2«l)2  I-F71  *n2**2«r?0*XL?**2: 

RUFI? I»G2*SQRT(G2**2-A. *r72*n?*Xl.?P|l 20TXl?-rZ|P02*F7?l  I  : 

DUFI  31*2.  *r.?-MIFI?l  ;  IF  ARSIBUFmi  <r.  ARSIRUM2H  TI1FN  IF  ARS  I IMIFI  ?1  1 
<  l.F-ll)  THEN  nil;  XI  3*1  .  ;H3«N2;GO  TO  SI  35FNO*.  n.SE  OOiO'IFIA  (  =  nui|?|  ; 

GO  TO  Sin;  F.NDS  ELSE  IHJF I  Al *RUFI  3 I  ; 

SIQ:  XL3*-2.*rZ2*D2/RUII4»:  l(3*XI  Ttll?: 

SIT!  7  3*H3f/ 2  S  IF  ARSIH3/73I  <  ACC  THEN  CO  TO  Sl9;  OISE  FZ0=F7.li 
F71 =F72 J 

•SI  60:  IllIF  111*1,1  F 22*  Al  NN*(  I T  PH  l-NN  VO  I  RY  -1;  RllFI  I  l*RHF  I  |  I  *7 3 5 
F22*F22«RI»ri1  l*AI  I  I  :  F.NO;  IF  ASSI  F7.2/F7  1  I  <  10.  TIIFNGO  TO  SIR; 

ELSE  XL3*XL3/?.1  H3=XI.3*N2;  73"H3«Z?;  GO  TO  S140: 

SIB:  »P*II35  XI 2*  XL3,  73=7(1  7  1*725  72*73:  D2=1.fXL2*.  GO  TO  S3; 

S19I  nP*.Mr»*l;  POOT  |i'R*=  7>:  IF  H*  *  N  THEN  00:  00  l«?  TO  NN; 

All  I  *  A I  1-1  I  v  i!  i« H  I N  r  I  •  A 1 1  I ;  END*,  M=r;N;  Gfl  HI  Ml  LNn;  ELSE! 

S25:  A*R(GIM;  IT  1CMK*0  TIIGN  RETURN;  FIST;  RFGEN=.D*A:  REGEN  I  11*1.1 
REGFN I  2  I  *- ROOT  I  111  00  K»2  TO  Nl  00  J=K*l  TO  2  RY  -11  PEGEMI Jl *REGFNI J» - 
ROOTI  111  ♦PEGFIII J-1  ll  run;  CNIlt  KfcGrN*AI  1»*HEGEN;  END  nCROliT  5 
11=11 

S=.01 5  PL  *100.;  OFT* 1 00. 5  RCT  :  Of)  KA*1  TO  10;  PL=PL  fDET*. 

no  JA*-S  TO  05 

R*  8 • AGE- 5*  S* JA*VS/77 .E ♦ 1 0  5 

Fl*-H*0.47*S**'5  F?=-|,  175E9*S**?5  T 3=M*1 . 5605F B*S**4 5 

F4»-l  ft.F 1 0* S*RFV.'*9.  CH>  *, 

F5«9  ,F.?n*l  |l,.’>SI**2-7.)4025F-9*S7P4»*N*2.7EIA*R*S-2.P2UE7: 

F6*-6. 34 7F-I  !i*S**2*PL  ;  F  7*  l-M*B . 607E-Bf2 .  B4C-B* S*R  I  *S**2*P L  5 
FR*I  M*3,  44?C3*R*S-1.904F6*SP*?1*S**2»P15  F9*-W»2.SB?F.  I  3*S*P4*PL  S 

R(1  1*1.5  RI2l*ri4T4:  Rl  3(*r?*Fl*r4FE0FF6; 

IH4I«F1*F5*F2»F4*F3fF71  M  51  *-F2*F5fF  3*F4fF  «H  Rl 61 «F3*F0«F9; 

NA=01  ICHK=0;CAU.  1>CP0(IT  I  0,  NA,  1  .  F- 1 0 ,  ROOT  ,  IClK,  RGI  5 

PUT  FILE  ISYSPRINTI  EOIT  |I'L|S,JA,RI  I SK  I  Pi  31 ,  E 1 1  1 , 41 ,  FI  10, 2 1  ,  FI  1 0 1  ,F  I  20,  1 0 1  I 
PUT  FILE  ISYSPRINTI  DATA  I ROOT  I ; 

nn  JJ*l  TO  55  IF  I  MAGI  ROOT  I JJU  <  0.  THEN  005  1F  11*3  THEN  110; 

PUT  FILr  ISYSPRINTI  FRIT  IPLI  IFIIOtAll;  GO  TO  LAI  TNI); 

ELSE  nn;  II=IIf11  PL  «PL-OET ;  DET=0»l*0ETi  GO  TO  RET;  ENO;  EN05 

else; end;  cno;  end;  la  i  eno; 


1.38 
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APPENDIX  B 


Difference  Operators 


For  the  computer  solution  of  the  laser-fluid  equations 
described  in  this  study,  the  differential  ecruations  were 
converted  to  difference  equations.  In  order  to  minimize  trunca¬ 
tion  errors,  very  accurate  seven  point  difference  formulae  were 
used  to  represent  the  differential  operators.  These  formulae  are 
somewhat  troublesome  to  construct,  but  a  useful  procedure  is 
described  below.  At  the  end  of  this  appendix,  a  few  of  the 
difference  formulae  used  in  the  computer  analysis  are  listed. 
Definition:  The  n  point  interpolation  polynomial  (abr.NIP)  of  the 
function  f(x)  in  terms  of  the  points  x^,x2«  . ..  ,xn,  is  denoted  by 
{f(x)}n  and  i»  a  polynomial  of  order  (n-1)  approximating  f(x)  in 
the  interval  [x^,xn].  {f(x)}n  has  the  property  that 

{f (x)}nlx-xis=f ^xi^  £or  i=1'2'  •••  'n  (Bl) 

and  is  explicitly  given  by  the  Lagrange  formula 


n-1^4  (x) 

{f(x)>  =  I  — rf  (x.) 

n  j-0? i(x.)  1 


(B2) 


where 

^(xj-fx— x^  (x-x2)  ...  (x— xi_1)  (x-xi+1)  ...  (x_xn)  (B3) 

In  other  words  «^(x)  is  the  product  of  all  factors  of  the  form 
(x-x^)  except  that  the  factor  (x-x^)  is  omitted.  For  convenience 
the  notations,  based  on  (Bl) , 


(f(x))„lx.Xi- 


{f  (xi))n“f  (xi)“fi 


(B4) 


will  also  be  used. 

Definition:  If  L  is  a  differential  operator  and  f(x)  is  a  suitably 
differentiable  function  of  x,  the  n  point  difference  cruotient 
(abr.NDQ)  representation  of  L,  denoted  by  is  defined  by 

the  following  relation: 

(  L)n  yf  (x)=L{f(x)}  nlx=y=L{f  (y)}  (B5) 

The  subscript  y  on  merely  indicates  the  point  at  which  the 

expression  is  to  be  evaluated.  As  in  the  case  of  the  NIP,  the  list 
of  the  n  mesh  points,  x^,  is  suppressed  in  the  notation.  For 
onvenience  the  notation 


where  k=l,2,  ...  ,n 


(B6) 


will  also  be  used. 

The  NDQ's  of  the  operators  bave  the  nice  property  that 


<f<x>}  «  "£l  iazYlifili,  }  f  (y)  (B7) 

n  J»0  J  i  dxJ  n,y 

for  any  x  and  y  in  the  interval  fxi/xnl»  since  the  Taylor  expansion 
of  any  polynomial  about  any  point  produces  that  same  polynomial. 

As  special  cases  of  (B7) ,  one  has 

f ,  -f  (x . )  ■"j1  jd_}j  f  (y)^!1  f(x)  (B8) 

1  1  j=0  j:  *dxJ  j=0  j.’  ‘dx  _  . 


n,y 


n  ,k 


for  any  i  and  k  in  the  set  1,2,  ...  ,n. 

Now,  by  combining  (B2)  and  (B5)  (or  by  inverting  (B8) ) ,  one  finds 
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<L}m 

ax' 


n,x 

n 

=  2 


•<§-> 

{f (x)}  «(S_ 

dx 

n  dx 

n. 

n,  ^(x) 

2  ** 
=1 

j  =1  W 

m 

m  n  V>.  (x) 

2  i - f. 


i  '"i 
f,- 


(B9) 


n 


where  the  prime  on  £  indicates  that  j,  ^i  and  j  ?*j.  ,  for  k^k# . 

-j  k  k  k 

k 

Now  (B9)  can  be  written  in  th«  form 


n 


{—-}  f(x)=  2  Aj(x)  f  *  where 

dx 


n,x 


i=l 


.  ,  .  n'  ^.(x) 

Ai  (x)  “  2  *  *  *  2  1 


j.  _i  j  1  *.  (x.)  (x-x.  )***(x-X.  ) 

^l"1  -*m  1  11  3l  3. 


(BIO) 


'm 


Thus,  if  one  agrees  to  use  the  NDQ  as  his  approximation  for  a 
differential  operator,  then  (BIO)  shows  how  to  express  the  m^ 
derivative  of  f(x)  at  the  arbitrary  point  x  in  the  interval 
[xl»XjJ  in  terms  of  the  values  of  f(x)  at  the  prescribed  mesh 
points  x^.  Actually  (BIO)  is  somewhat  more  general  than  required 
for  the  computer  program  described  in  the  paper.  For  the  purposes 
of  this  program,  only  the  following  special  case  of  (BIO)  is  reauired: 


its £*-  * 

' ax,n,k  i=l 

*i  k"'  s’  7—^*7 - T-S-; - r  (BID 

1,R  n_1  3m_l  *i(xi>  <xk-xj  >"’(xk-xj  > 

-*  1  -'in 


where  i,k=l,2,  ...  ,n  and 

_  1 
I  • 

’1  "m 

Formula  (BIO)  or  (Bll)  guarantees  explicit  NDQ  representations 
of  any  differential  operator  L  and  shows  a  way  to  compute  it. 
Although  a  great  deal  of  computation  remains,  even  with  formula 
(Bll),  at  least  the  work  is  straightforward.  Clearly,  the 
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definitions  and  procedure  shown  here  are  trivially  generalized 
to  cases  with  more  than  one  independent  variable. 


In  the  list  below,  the  various  difference  operators  used 
in  the  laser-fluid  computer  program  are  given.  For  each  indepen¬ 
dent  variable  the  grid  spacings  were  taken  to  be  uniform  with 
sizes  denoted  by  h  or  k.  The  great  variety  of  formulae  was  used 
in  order  to  handle  the  special  situations  near  boundaries. 

A.  Center  Difference  Scheme 


(B12) 


T 


{axK»4fiSh2[2*fi+i+fi-l)  20(fi+2+fi-2)+50(fi+3+fi-3)  18fi] 


B.  Off  Center  Difference  Scheme 


1. 


{^}7,4+l£i=1h['T2£i+I£i±l'ifi:':l‘,s£i±2'kfefiT2*l5fi±3~|o£i±4] 


upper  signs 


*  &  * 


lower  signs 


[  10fi4rlfi±l' 

”iofi±5] 


107, 


17, 


11 


"90ri+l‘tT6ri±2"l80fiT2  10Ii±3 


rf. 


45  i±4 


(B15) 


upper  signs 


lower  signs 
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(B13) 


(B14) 


2. 


\  Ox 


"}7,4+2fi=i^["  Jsfi+  Ifi±l"  i±2+  ffi±3"  ifi±4+  Iofi±5] 


1 

upper  signs 


lower  signs 


ij_\2  f  J  L  7  f  .  27f  .  7  f  .  19f  _  67f 

Ox/7'4*2  *  h2[  18  1  10  1*1  10Ein  4  i*2  18  i±3 


+  |f .  ,  A 

5  i±4 


-If  +  11  f 
2  i*5  180  i±6 


j 


(B17) 


10- 


upper  signs 


C.  Mixed  Differentiation  (Center) 


lower  signs 


{dxdy}fij"2Ek[2(fi+l,  j+l+fi-l, j-l"fi+l,j"fi-l, j"fi, j+l”fi,  j-1* 

“  2^(fi+2, j+2+fi-2, j-2“fi+2, j“fi-2, j”fi, j+2“fi, j-2>  + 

4  90(fi+3, j+3+fi-3, j-3“fi+3, j”fi-3, j”fi, j+3”fi, j-3* 

+  il£ij]  I  ■  <B18> 


D.  Mixed  Differentiation  (Off  Center) 


I  a2  If _ l_r  21f  ,  13r  .  107-  .  17- 

Idxdy/  13"  2hk["  10Eij+  I8£i+l,j*l+  90£iTl,j±l+  36ti±2,jq:2 

_  11  f  _  _  3 f  _  .  4 f  _  1 f 

180  i^2 »  j±2  i±3,  j^3  j^4  gg  ii5»  j^5 

-h2(^?fij-k2(^,2£iJ  (B19) 


(B16) 
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17  27  7  19 

dxdyi  ij  2hk  18^i  j  10^i±i#j^i+  10^i^i,j±i+  4^i+2,j:f2 


-  &k. 


18fi+3,j*3+  5fi+4,j+4"  2fil5,j+5+  180fi±6,j+6 

-h2<|j>2firk2(^,2£i:] 


(B20) 


•  i 

•  i 


upper  signs  lower  signs 


(Additional  points  on  the  axes  are  to  be  selected  for  the 


unmixed  derivatives) 
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APPENDIX  C 


COMPUTER  PROGRAM  FOR  THE  LASER-FLUID  EQUATIONS 


In  this  appendix,  the  computer  program  for  the  laser-fluid  system 
used  to  obtain  the  results  discussed  In  Section  IV  Is  exhibited.  Before 
programming,  certain  simplifications  were  made  In  Equations  (174)- (186) • 
based  on  the  fact  that  a,  the  linear  absorption  coefficient.  Is  very 
small  In  air  at  10.6  u.  For  x  f  0,  Equations  ( 1 74)- (179)  were  reduced  to: 


l£  =  .  Jp[ib.x)2  +  (1-x),  v  +  +  XI 

3t  r0  3x  rQx  r  3z  J  r 


~x)  y  +  y 

r_  r  3x  z  3zi 


(Cl) 


3T  _  1  T(l-x)4  32T  .  (1-x)3  / 1  -2x\  3T  .  32Tl 

3t  P  ( c;  L7^"  T?  7T”  [—) 


+  2n+nl 

Lv 


ro  x 


r 


,  T/i  i3  3V  .  3V  #,  \2  3V  3V  1 

Hi  (V-x).,,  y  _ r  +  ]^x  y  _z  +  (1-x)  _V  _lz 

Cv  L  r  *  x  r  3x  rox  r  32  ro  9x  32  J 


v  l  rQ  x 


*  'i'1-”  [S  *  ‘Jf1-  v  If  *  <■  8]  •  %  *  1  '*’) 

•f 


o-xr  y  ii+v  ii 
r  r  3x  z  32 1 
0 


(C2) 
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For  x 


3Er 


8t  e 


1  t  c2 

r  ^ 


(1-«)4  til  I  (l-«)3  (hi*)  8E1  t  !lii  +  ;k, 

.  r  ax^  r  Z  '  x  '  3x  azZ  L  3z  j 


[s[E"(!t)  *t  (ve)]  Ei  •  c‘  It  e2} 


0,  equations  (181)  -  (186)  become: 


H  If  ] 

0 


avr  av. 


r  r 
o  o 


+  ^-^(e2  +  e2)J-vz|1 


(C6) 


(C7) 


(C8) 


3V 

t 

3t 


(C9) 


lj-§TH+(W) 


+ 


(CIO) 


(CD) 


1 

4 

(Cl  2) 

In  this  computation  the  following  notations  have  been  Introduced  Into 
the  computer  program: 


N 


The  number  of  steps  In  x  (or  r)  to  be  remembered. 


NN 

• 

• 

Total  number  of  steps  In  x. 

M 

• 

• 

The  number  of  steps 

In 

• 

MT 

• 

• 

The  number  of  steps 

In  time. 

DX 

• 

• 

Ax 

EZ 

• 

• 

Eo 

WP 

3<V1) 

DY 

• 

• 

AZ 

KL 

• 

• 

*1 

WQ 

_l 

csj 

o 

DIT 

• 

• 

At 

KM 

• 

• 

2KL 

WR 

c  /2u>l 

VI 

• 

• 

X 

WA 

• 

• 

2nV 

WS 

1/2u)l 

V2 

• 

• 

<l-x)/r0 

WB 

• 

• 

nW 

WT 

Rp/M 

V3 

• 

• 

(l-x)2/r0 

wc 

• 

• 

n 

WU 

k/Cv 

V4 

• 

• 

(l-x)3/r02 

WO 

• 

• 

n' 

WV 

2 

V5 

• 

■ 

<l-x)4/r02 

WE 

• 

• 

t(y-i) 

WW 

(2n+n‘) 

V6 

• 

• 

1/x 

SF 

• 

• 

-RT/M 

WX 

ac/2Cy 

V7 

• 

• 

1/x2 

WG 

• 

• 

1/ro2 

WY 

1/2 

V8 

• 

• 

(l-x)/r0x 

SH 

• 

• 

2/r-o2 

WZ 

^o 

V9 

• 

• 

(l-x)3(l-2x)/r2x 

WI 

• 

• 

po 

VA 

2/r0 

VI 0 

• 

• 

(l-x)2/r02x2 

WJ 

• 

• 

eo 

VH 

n/Cv 

DU 

• 

• 

1/Ax 

WK 

• 

• 

V1 

VI 

2n'/Cv 

DV 

1/(AX)2 

WL 

• 

• 

Wj_/2 

VN 

1 
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OW  :  1/Az 

*  ! 

VO  :  273 

DZ  :  l/(Az)2 

MN  :  3(e0+2)(£0-l) 

VP  :  3 

DT  :  l/2(Ax)(Az) 

MO  :  2n/r02 

VE  :  (y-1) 

HI  :  l/C(e0+2)p0- 

H3  :  e 

H6  :  1/p 

(e0-l)p] 

H4  :  | 

H7  :  1/e 

H2  :  (e0-l)  HI 

H5  : 

dp2 

H8  :  /e 

H9  :  (eQ-e) 

X(l.K.L):  E, 

X(2,K,L):  E? 

X(3,K,L):  p 

X(4,K,L):  T 

X(5,K,L):  Vr 

X(6,K,L):  Vz 

PXA.FXB.FXC:  |j- 

32 

SXA.SXB.SXC: 

3  X* 

SZA,SZB,SZC: 

32 

FZA.FZB.FZC: 

32 

MDA.MDB.MDC: 

Following  Is  the  computer  program  for  this  study. 


t NOUNOERF LOH )  t  OPTIKSl  PROCEDURE  OPTIONS  (MAIN)  t 

NN-23I  M«8lt  MT-lUOt  N-16;  MTT-lOOt  ZNN-NNt  ZN-Nt  ZKT-MT I 

OCL  (LY,S1<7),S2<7),S3(7),S4(  7I.S5I7I ,S6«7)I  LABEL! 

OCL  CCZ.UN,  VLtVr,  CMC.  ZIJ.ZZ. 

UX.UU.OV.OY.UW.UZ.UT.OIT  (OA .Ob. OC .00.01:. OF . UC.OH.Ot  ,  (iJ.OL  ,OM,ON.OO,OP.  00, 
0K,UP,H8,HQ,Hl3,hll ,H12»Hl3,Hl4,PHT,Gl( 121) t 
XR.  O&.DT • OUt  UV. OH.  OX.OV .DZ.PA .PM.PC .Pli.Pfc.PF . PG.PM.PN. 

HA. Hll  . HC «  HII.HC.HF  .ViO.WII.HI  .HJ  .hK.HL  .hm.hn.hu.hp.hu. hh. 

HS.HI  .HU, WV.HH.HX , HY, HZ .HI  .HZ  «H  3,  W4,  H5,  H6  «t<  7 ,  VA.HA.HB,  III).  ME ,  HF ,  ItH  ,HS,H  1 ,  M2.  HJ, 
H4,H5,lt7,HB,H9, VH, VI ,VN,V0,VP,V1,V2, V3,V4,V5,V6,V7, VU,VV,V10,F12, 
KL.KM.Xl.X2.X3.X4.XS.X6.X7.XU.XV.X10.F2,F3.F4,FS,F6tF7tF6,F9,F10.Fll.A(6).B(6). 
C<6), 0(6), 1(6), M6),X<6.16, 121), Y(6,  16,121))  BINARY  FLOAT  (  53)  i 
OCL  (FXC.SXC.FXA.SXA.FXU. SXB.MOC .MOA .MOUeFF ,F FI )  RETURNS  (O)NAKY  FLOAT  153) )t 
EP-1.L3;  I J* ICO!  Z1J«)JI  CHG-l.I  IC-l! 

NA-N-lI  Nil-N-2;  NC-N-31  N0-N-4I  N5-N-5;  NF-N-6I  N7-N-7I  NNA-NN-1 I 

MA-M-1;  MB-M-2)  MC-M-3I  M0-M-4I  ME-H-5;  liF-M-61  HG«M-7t  ZMA-MAI 

UU°NNA!  UX-l./MUJ  OV-l>U*OUI  ()Y  =  3 .  6C  6/  ZHA I  UW-Z MA/3. 6E6 I  OZ-OH*OH; 
0T=.5*UU»DHt  01T-1.L-.5/ZMT!  0P-1./1H.;  Kl  =  SUKT  (  1 .  +5. 65 1-4)  *5 . 91 E3J  KM-2.*KLt 

UA-.7SI  III)— .151  DC  *  1 .  /  6C .  I  OU-1.5!  OE-1./90.J  UF— 4V./18.I  0G-2.lt 

011*13. /18. I  01*107. /VO. I  UJ-17./36.I  UK— 1.1*0P;  01.*-. 3l  0M-4./45.I 

ON— 3LI  00--7./1B.I  0P»-2.7i  0Q-.7I  XR-4.75;  US— 67. /IB.  I  OU— .51 

OV--UK;  OH-l. 5*001  (IT-1. 01  OX-4. /3. 1  0Y--.4,-  UZ-2.*UCt 

PA*1.5*UMI  PC— 77./60.I  PO-2.5I  PE  —  1./6.I  PF— 5./3.  I 

pc-5.*rc;  Pii- .251  Pt;— OX'-OXj  J’M—  OYFCYi 

HA-3.55E-4I  HO* 1 • 775E-4  I  HC-HBl  HV-2.I 
VP-3. 1  HC-O.t  VC*. 4  I  HI-1.25E-3S  HJ*1 .♦5.65E-4 1  HK-5.65E-4S 

HZ-5.L-3I  HG*HZ*WZ I  MH-2.-WGI  WL-.UB65LI4;  VF— 8. 3144E 7/29.01 2; 
HM-I  HJ*HV)  *WI  I  HN-VI'*HK*WKI  HO-Wt*wHt  M**1.69SE-3|  W0*9 .Eo/ 1 . 773  I 

HR*. 5*H0l  HS-.5/ 1. 773C14I  HT— VFl  HU*. 1B6*H) I 

HK-WA/ 7. 1431:61  HX-V.E-3/1. 42061  KY-,51  VA-2.*HZ) 

Vll-HC/  7. 143C6  •  VI  *0.1  VN-I.S  VO-283. 1 

HI— OV*  134 . 89  / 1  8. I  H2-0V*12.  I  W3— 0V*7.5t  H4»0V*40./V. I  H5--0V* 1. 6 75 1 

H6-0V*.48I  h7— OV*OPI 

H6*0U/(441.+1RC.*UX)I  HQ-1080.  *hli  i  H1C— 675.*H0S  HU-400. *H8t 

H12  — 160.75*N8t  H13-43.2*H0I  W14  — 5.*HUt  H8— 674 ,4b*HU  I  RT I 

EZ-(2.C2/Cl.*2.82E-4))**.5/3*EP5 

00  K-l  TO  N;  00  L“ 1  TO  Ml  ZK-K;  ZL“L{  UN- I ZK-VN ) /OUt 
X(l,K,L)=CZ*EXP(-(rtV*UN/( l.-UN) )**2-( . 10*ZL-3.0)**2) l 
X(2,K.L)*0.t  X(3,K,L)*H);  X(4,K,L)*V0S  X(5,K,L)*0.t  X(6,K,L)=0.i 

ENOt  LNOI 

00  L-H  TO  1  UY  -It  PUT  FILE  (SYSPKINT)  EUIT  (L, IK, XI I ,K.L)  UO  K*1  TO  NA)) 

I  SKIP (2). XI 4). FI  3)  ,4(SK)P,5(X(5.' ,F(2),X(2),E( 14,7))));  ENUI 
FXCS  PPOC  IX1.X2.X3.X4.X5.X6)  01 NARY  FLOAT  153)1 

UCL  IXl,X2,X3,X4,X5.Xh)  ID  NARY  FLOAT  153)1 

F2»OA*|Xl-X2)*O0*(X3-X4)+OC*IX5-X6)l  RETURN  (F2>;  ENU  FXCJ 
SXCS  PRUC  (Xl,X2,X3,X4,X5,X6.X7t  BINARY  FLOAT  153)1 
OCL  (X1,X2,X3,X4,X5.X6.X7)  II I  NARY  FLUAT  (53)1 

F3-00*(X2  +  X3)*ab*CX44-X5)4-l)E*(  X6*X7  )»OF*XU  RETURN  IF3)I  END  SXCt 
FXAI  PROL  (X1.X2.X3.X4.X5.X6.X7)  BINARY  FLOAT  (53)1  OCL  ( X 1 , X2, X3 , X4, X5 , X6. X7 ) 
BINARY  FLOAT  (53)1  F4 -|'C*X1 »P0“X2 «PE*X3U'F*X4*PG*X5*PH*X6*UZ*X7 ;  RETURN  (F4II 
ENO  FXAI 

SXAl  PRUC  (XI ,X2,X3,X4,X5.X6,X7,X8)  BINARY  FLOAT  (53)1  OCL  I XI, X2, X3, X4 , X5, X6 , 
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X7 1  XB )  It) NARY  FLOAT  (5311  F5«UO*Xl ♦UI»*X2*O0*X3*XR*X4*US*X5+OT*X6*0U*X7*UV*X0) 
RL1URN  0  5 )  |  END  SXA; 

FXUi  PROL  (XI  ,X2,X3,X4,X5,X6, X7I  BINARY  rLUAT  (5315  DCL  ( XI ,X2, X3 , X4, X5, X6, X7I 
BINARY  FLOAT  (531)  r6«0W* XI ♦0X*X2*0Y* X3*UU* X4*07*X5*PA*X6-0C*X7 5  lit TURN  (F6)) 
END  (  XII 5 

SXUI  CROC  (XlfX2fX3iXAfX5«X6«X7|Xd)  BINARY  FLOAT  (5315  OCl,  ( XI , X2 , X3, X4, X5, X6 , 
X7.XIU  BINARY  FLIIAT  (53)5  F7*0G*X1  ♦UH*X2*(U  *X3*0J*X4<  OK*X5«-UL*X6*UM*X7*UN*XOS 
RETURN  ( F 7 ) 5  END  SXU5 

HUCI  )>IU)C  (Xl«X2iX3tX4tX5.XbiX7«XH|X9)  BINARY  FLOAT  (53)5 
OCL  (XI tX?tX3«X4fX5tXb,X7tXBlX9)  BINARY  FLOAT  (53)5 

F8*0F  *  XI  ♦)!))*)  X2*X3)*l)0*(X4*X5l*O£*(X6iX7)*PH*XB*PN*X9s  RETURN  (F(l)5  END  HOC  5 
HUAI  )’K(IL  (XI  lX2tX3fX4,X5«X6lX7f  XU|X9IX1C)  (UNARY  FLOAT  (53)5  OCL  (  XI ,  X2 .  X3  ,  X4 , 
X5,X6,X7,Xil,X9f  XIO)  (UNARY  FLUAT  (53)5  F9*U0*Xl ♦OP*X2*GQ*X3*XR*X<UUS*X5«OT*X6* 
0U*X7*(iV<‘XII*l*M*X9*PN*X10S  RETURN  (F9)5  ENO  HIM) 

MORI  )>R()C  (  XI  (X2iX3tX4t  X5iX0,  X/ ,  XB ,  XS  ,  X I  0  I  BINARY  FLOAT  >531!  OCL  (  XI  t  X2f  X3 ,  X4, 
X5,Xl»iX7,xelX9,X10)  Ml  MARY  FLOAT  (53)5  F  10=()G*X1  ♦0H*X2  *U1  *X3+0  J*X4*0K*X5*0L*X6» 
UM*X 7 ♦UN*XMO,M* X9O’N-»X10 5  RETURN  (F10))  ENO  HOBS 
IF:  CROC  (XlfX2,X3.X4tX5«X6tX7)  BINARY  FLOAT  (53)5 
OCL  (XltX2«X3lX/ifX5fX6fX7l  BINARY  F'LUAT  (53)5 

Fll=>W)*XUw2*X2»H3«X3*'.*'.'eX'.*W5*XS*W6*Xb*WOX7)  RETURN  (Fll)!  ENDFF! 

Fl:l  8  I'KOC  (XI  ,X2,X3iX4,X5,X6,  X7)  BINARY  FLOAT  (53)5 
)l(.L  (  X)  ,X2,X3,X4,XS,X6,X7)  II )  NARY  FLUAT  (53)5 

F12'rfR*Xl4W9«X?  +  WlC*X34Wll<'X<i«nl2*X5  +  WlJ*XbUm*X7SKCTURN  (F12I5  ENO  FF1S 
OU  10=1  TO  HT T) ) F  )0  >  1  THEN  U0> 

IF  1(»IJ  THEN  00;01T=5.E-0;  Z7  =  3.L3*(/.IJ-1.I*DW!ZT*1D-1J*1S  CHG=0.5!  fc NO 5 

ELSE  UN 5  /7*0. 1  71*105  ENU! 

00  K*1  TO  N!  ZK«KS  ON*)7K-VN)/DU5 

00  L* 1 1  M)  ZL*L8 

X(l,X,L)-E2*EXIM-(WV*0N/(  l.-ONI ) ♦»2-( . 10* ( 7L-30.-77-3. E3* ( 7T- 1 . ) ♦CHG*DWI )**2I ) 
ENOS  END) 

DO  1*2  TO  MA!  2L*L  5  ON* ( 7N- VN ) /OU! 

X(liNtLI*F  7  *F  XI’  (  - )  RV*ON/(  1.  -UN )  )**2-(  .  10*(ZL-30. -Z/-3.F3*  ( 2T-L  .  )  *CIIG*I)W|  )*#2 1 1 
ENO!  tflOS  ELSE) 

0(1  L=  2  TO  M- 1  5  00  K»1  TO  N-l 5  ZK=K1  00  J*1  TO  6!  A ( J I  =  X( J , K,L 1 5  END! 

II7*WI’*A(3)  5  HI*  VII/ 1  WM-KK*A(  3 1 1  5  H2  =  WK*lll!  H3*VN»H7*Hl5  H4  =  WN*IU*Hl! 
II5*WV*II4*H?  5  H7*Vf,-M7/)WM*|,  130E-3*A(3I )! 

110*1)3* *0.  5!  H9=WK*(3.«hK)*)»il-AI3l  l*Ml  J 

IF  K  >  3  THEN  00! 

V) * ( 7R- VH ) /OU!  V2*)VN-V1 )*W7 5  V3  =  V2*( VN-V1 1 5  V4=V3*V2;  V5*V3*V3! 

V6*liU/ )  7K-VN I  5  V/=>V6*V6SVU»W2*(2NN-ZK)/(  74- VN)  5  V')*  V4*  ( V6-WV )  5  V1C*VB*V85 

IF  K  <  NC  THEN  00!  KA*KMS  KB*K-1)  KC*K«2! 

KO*K“ 7 5 Kl. 3 5  KF*K-3!  DO  J*1  10  6! 

))(J)*I:XC)X(J,KA,L),X(.I,KB,))#X(J,KC,L),X)  J,KO,L),X(J,KE,L1,X(J,KI,L)I*OU!  . 

C(  J)*SXt) A( J) ,X(J,KA,L) tX(J,KO»L) ,X( J.KC,LI,X( J,XO,l  ),X( J,KE,L>  , X ( J ,KF , L ) I *0V! 

ENO!  IF  K  >  7  THEN  LX»7)  ELSE  LX*65  END! 

ELSE  OU!  >1  K*NA  IIDN  )>(IS).X*58  UO  J=1  (0  (.!  B ( Jl  =-FX A)  A (  J ) , X) .1 , NO  ,L  ) , X(  J , N, L ) , 
X(JtNCiL) ,  X)  JtNIlf  I.)  tX)J,N‘.,L)  ,X)Jt  HI •.(.  I  )*.)U!  C)  J  I*  SXA )  A  (  J I  ,  X)  J  ,NU  ,L ) ,  X(  J,  Nt  L  I , 
X  ( J  (NC  »t  )  t  X (  J  tNOtL  )  t  X )  J ,N5*(. )  f  X)  J ,NF# L  )  *  X(  J*N7,L I  l*UV;  END!  END;  ELSE  DO! 

LX*  4!  DO  J*l  10  (.|  U).l)=-FXB)A)J),X)J,NCtL)iX(J,NA,L),X(JtNO,L)  *  X  )  J  *  N .  (  ). 

X(  J  *N5»L  )  t  X  (  J» Nl  •  L  1 1  *0U !  C 1  J 1  =  SXB>  A 1 J 1 ,  X(  J. I1C  i L  1 ,  X ( J ,  NA ,L  1  ,  X(  J , NO  ,L  1  ,X(  J ,  II, 1 1 
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f  X(  J  , NS  «L  I  ,  X ( J , NF , L I  « X ( J , H7 , L I  I *DV;  ENOS  END;  END;  END;  ELSE  DO; 

IF  K»1  THEN  DO;  LX«lJI)l)  J*1  TO  4  , fa ;  It  (  J  1  “0 .  !  C(  J  I  bFF  ( A  (  J  |  «  X ( J  ,  2  ,  LI  ,  X  ( J  ,  3,  L  I « 

X(  J,4,L)  ,X(  J,5,LI  ,X< J,fa,l I  , X(  J ,  7 , 1. 1 1 ;  END; 

(Mb)  »Fri(  A(S)  ,X 1 5, 2,1 1,  X(!>,3,L),X(5,4,L),X(5,5,l),X(5,b,Ll,X(5,7,Lll; 

C  (SI»2.*U(;>I  ;  ENOS  ELSE  i)U.‘ 

VI  * ( ZK-VNI /DU;  V2«(VN-VI l*WZ;  V3=V2*( VN-V1);  V4*V3*V25  VS=V3*V3; 

V6*0U/(ZK-VNI S  V7*Vfa*V6;Vt!“WZ*(ZNN-ZK)/(ZK-VNI  !  V'.»  =  V4*(  Vfa-NVI  !  V1GBV8*V85 
IF  K*2  THEN  DO)  LX*2;  00  J=1  TO  6  ;  tfl  J  l=FXA(  At ) ,  X  (  J,  3  ,L  I  ,X(  J  .  I  ,  L  »  ,X  (  J  i  4  i  LI  i 
X(J,S,LI ,X( J,b,L»,X(J,7,L )I*DU;  C( J 1  =  SXA< A( Jl ,  X ( J , 3,L) • XC J , l , L»  ,  X ( J  ,4, LI t 

X(J,0,LI  ,X(J,fa,LI  ,X«J,7,L).X(  J,8,L)  l*DV;  END;  END;  ELSE  DO.*  LX*3; 

Oli  J«>  TO  o;  (t(J)BIXG(A(J),X(J,4,L),X(J,2,L),X(J,5,L),X(J,l,L),X(J,6tL), 

XI  J,  7,  L)  l*UU;  C(  JI*SXIMA(  Jl  ,X(  J ,  4 ,  L 1 1 X  ( J,2,LI  ,X  C.J  ,S,ll  ,  X I J ,  1  ,Ll  ,X  ( J  ,6,L  I  • 

X ( J  *  7  T L )  ,  X(  J  1 0 * L I  )  *> D V ;  LNU;  END;  END}  END;  ‘ 

)F  L  >  3  THEN  DO;  IF  L  <*  MC  THEN 

oo;  la*l*i;  L(i*L-l.;  lc=l«2s  L0*l-2;  i  o*l^3;  lf«l-3;  oo  j*i  to  6; 

0(J)*FXC(X(J, K.LAI ,X(J,K,LOI,X( J.K.LC I ,X(J,K,LDI,X(J,K,L5I ;X(J,K,LF  1 1 *DW, 
t(JI« 

SXCI A( J) ,X( J,K,LAI ,X(J,K,U)I  ,X( J,K,LCI tX l J, K,LD) ,  X(  J,K,L5I ,X ( J , K , LF 1 1 *DZ ;  ENDS 
IF  L  >  7  THEN  DU! 

GO  TO  S1(LXI  ;  Sim:  LY-SA;  GO  TO  SK;  SII2):  LY»SS4;  GO  TO 
SR;  SI (31 :  LY* SS5 ;  GO  10  SK;  S1(4|:  LY-SS/;  GO  TO  SK;  SMbi:  LY*SS6;G0  TO  SK; 
sl(6l:  LY«SSl;  GO  TO  SR;  SI ( 7 1 1  LY*$S1;  SIC  ENDS  ELSE  00;  GO  TO  S2(LX>; 

S2(  1 1  :  LY«SA;  GO  TO  SO;  S2I2I:  LY-SS2;  1.0  TO  SO;  S2(3|:  LY=SS3J  GO  TO  SO! 
$2(41  :  LY«=SS9,*  GO  TO  so;  S2(5H  LY=SSO;  GO  TO  SQ;  S2(6) :  I.Y^SSU  GO  TO  SO; 
S2(7l:  LY-SS1;  so;  END;  end;  llsl  do;  ie  L“iiA  then  do,  do  j-*i  to  6; 

0(JIB-FXA(A(JI,X(J,K,KI1|,X(  J,K,MI  ,X(  J,*,HCl ,XIJ,K,MDI  ,X! J,K,ME)  ,x( j,k,mfi i*dh; 

E(  J>»SXA(A(  JI,X(J,K,HU).X(J,K.H)  ,  X(  .*  ,K  ,MC  I  ,  XI  J  ,K  ,MDI  ,X  ( J  ,K  ,KEI  ,X(  J  ,K  ,MF  |  , 
X(J,K,MGI l*DZS 

END;  go  TO  S3 (L XI ;  S3(ll:  LY-SA;  GO  10  so;  S3(2) :  LY-SS4;  GO  TO  so;  S3(3): 

L Y  =  SS 4 ;  GO  TO  SU;  S3(4):  LY»SSb{  GO  TO  SO;  S3(bl:  LYBSSbi  GO  TO  SO;  S3(6I: 

L Y*S54 ;  GO  TO  SU;  S3 ( 7 1 X  LY»SS6!  SO:  END;  ELSE  DO;  DO  J»1  TO  6; 

0( JI«-FX«IA( JI,X(J,K,HCI ,  X(J,K«MA),X(J,K«HD),X(J,K,M),X(J,K,HC I ,X( J,K, HEI )*DW; 

E( JI*SXO(AI JI,X( J.K.HCI,  X (J.K.MAI.XIJ.K.HDI .X( J,K,KI ,X^  J.K.Mfc) ,X( J.K.MFI , 

X( J,K*NGM 'DZ {  END;  GO  TO  S4(LXl;  34(11  :  LY=SA;  GO  TO  SH;  S4(2»:  LY*SS4; 

GO  TO  SH; $4(3 1 :  LY»S$5;  GO  TO  SH;  S4(4):  LY-SS7;  GO  TO  SM;S4(b):  LY«SSfa!  GO 
TO  SH;  S4 (6)  :  LY*$S5i  GO  TO  SM;  S4(7|:  LY«SS7;  SH:  END;  ENO;  END;  ELSE  DU; 

I F‘  L-2  THEN  00;  DU  J*1  TO  6; 

0(J1=FXA(A( Jl  ,X(J,K,3I,X( J,K, 1I,X(J,K,4I ,  X(  J  ,K ,  51 ,  X(  J ,K ,  6 )  ,  X( J,H,7I l*OW; 

E (Jl *SXA( A( Jl ,X(J,K,3lfX(J,K,ll,X(J,K,4),X(J,K,SI,X(J,K,6),X(J,K,71,X(JtK,8)|*U 
l  ;  END;  GO  TO  S5ILXI;  S5I1I:  LY*SA;  GO  TO  SK;  S!>(2|:  LY=SS2;  GO  TO  SK;  SSI3I: 
LY*SS2;  GO  T(l  SK;  SS(4»:  LY=SS8;  GU  TO  SK;  Sb(!>|:  LY*SS8;  GO  TO  SK;  S5(6|: 

L Yb SS2 ;  GO  TO  SK;  S5(7|:  LY-SS8;  SK:  END;  ELSE  DO;  00  J*l  TO  fa; 

0(JI-FXU(  A(  JI,X(J,K,4»  ,X(  J,K,2I  ,X  (  J ,  K  ,!>  I ,  X  ( J  ,K » 1 1  ,  X  ( J  ,K  .61  ,  X  (  J ,  K,7 1  l*OW; 
E(J)*SXO(A(JI ,X(J,K,4),X(J,K,2I .X ( J , K ,5  I , X( J ,K , 1 1 , X( J ,K ,6 1 , X ( J , K,  7 1  ,X( J  ,K  ,81 1 *0 
l  ;  END;  go  to  S6(LXl; 

S 6 ( 1 1  :  LY-SA;  GO  TO  SI;  S6(2II  LY=SS2;  GO  TO  SI;  S6(3|:  LY-SS3;  GO  TO  Si; 
S6(4I  i  LY-SS9;  GO  TO  SI;  Sfa(bl :  LY^SSO;  GO  TO  SI;  S6I6I  :  LY  =  S  S3 1  GO  TO  SI; 

S6(  71  <  LY=SS9S  Si:  END,*  ENOS  GO  TO  LY;SS1:  1)0  J»5  TO  6; 

F(JlcKUC(A(JI,X(J,KA,LAI,X(J,KB,LUI ,X(J,KC,LC),X(J,KD,LOI , X ( J , KE , L5 I , X( J,KF ,LF I 
,C(  Jl • E ( J I l*DT;  end;  GO  TO  SO;  SS2:  DO  J«5  TO  6; 


F  (J  )  =  HOAI  A(  J)  ,X(J,K«1,L*1  I.XI  J.K-l.L-l I tXI J,K»2,L+2)  ,  XI J fK»3 , L  +  3)  , X I J , K*4 ,L+4 1 , 
X(  J  fK»btL»M  ,  K<  J,  K«6,L+6)  ,C(J  )  ,t(  JI1*0T  {  END;  GO  TO  Stii  SS3:  00  J=b  TO  65 
F(  JI=MOBIAI  J)  ,XIJ,KtliLMI,X(  J ,K- 1 , L- 1 ) , XI J , K*2 ,L+2 I  ,X(  J,  K-2.L-21  ,X  I  JfK< 3,L*3) t 
X(JiK*A,L  +  6)«X(J,K<5«L*5)ft(J)it(JH*OT;  END;  GO  TO  SO;  SS4:  00  J  =  5  TO  65 
F(  J1=-MDAI  A( J) ,X( J,M l,L-H  tX( J.K-l.LU 1 . X( J , K»2, L-21 • X( J.K ♦ J.L-3J * X( J ,K»4,L-4l 
» 

Xf  J  »K»E-  ,L-5>»  iXI  J,M6,L-6I  ,CIJ1  ,f.U|  l*OT;  CKO!  GO  TO  SO;  SSS:  00  J  =  S  TO  6; 

F« J)=-HLI0«A« J),X(J,K*1,L-1),X(J,K.  l,L*I  1 1 XI J ,K*2, L-2 1 , X I J , K-2 , L+2 1 , X( J , K  +  3, L-3) 

i 

X«  JiK+4,L-4)  t  X(  J,K*S,L-5t  »CIJ1  tEI  J1  l*DT;  ENO;  GO  TO  SB;  SS6:  DO  J  =  5  TO  6} 

FIJI  “HO  A  (  A  (  J 1  ,XIJ,K-l,L-ll,XIJ,K+l,L*ll,X<J.K-2,L-2»,XIJ,K-3,L-3)tXIJ,K-4iL-4), 
X(JiK-!»,L-bl,XIJ,K-6,L-6T,CIJ),EtJTK0T;  LOO;  GO  TO  SB;  SS7:  DO  J  =  5  TO  6; 

FIJI  ■MODI  A( J I iXIJ  rK-ItL-i I. XI J  ,K+ i »L  + 1)  .XI J,K-2,L-2l ,X I J ,K+2, L>2 1  ,XI J , K-3 ,L-3)  , 
XI  J,K-4iL-4)-,XI  JiK-5,L-!i),t(J)-tl(J't)*OT5  EMI) ;  GO  TO  SB-;  SS8 :  00  J  =  b  TO  6; 

FI  JI=-MOAI AIJ l,XIJ,K-l,L+l) ,XIJ|K+1,L-1 1 , X  I J ,K-2 , L« 2 > ,  XI J , K-3.L  +  3) ,XI J,K-4,-L+41 
» 

•X  I J  »K-b»L*  b  I  «  XI  J  «  K-6,  (. «  67  it  I  J  1  i  E  I  J  1  1  *UT  S  ENDS  GO  *T0  SB;  SS9»  DO  J=*5  TO  '6  5 

I- 1  Jl  =-M»»»l  Al  J  I ,  X I J  ,  K-  1 ,  L  ♦  II  ,X(J,K+l,L-l),X-<J,K-2,L  +  2l,X<J,K<-2,L-2liX<  J,K=3,L*3) 

XI.J.,K-4,L*4I  ,XlJiK-5t(  +51  ,C(.)  I ,  K  J  1 1  *01-5  UMt();  GO  TO  SB-; 

'SAS  IIA-AI  6  I  #0  I  3  I  ;  1111=  AT  IT  #AI  1»H  Al  21*AI*l!5:  VI  3t  K-,  1. 1 «- I HA+A  I  31  *  I  D(  6  T^VA^B  1 5)Vl  5 
IIO=ll,.>*YI3.K,LH‘YI3,K,LT  5  III  =  H4*Y'(  3  ,  K  ,  LVi  'HH=WS*HU+WL*H9  5 
Mt  =  Vl;  *  A 1 4  I  5  KK=VI:*Al4l5 

Y  14  •  K  i  L  )  ■  I  rlU*  I  Wit*  C  14  I  tfc  1 4  1 1  ♦MW*I  01  6 1  *01 6  I  ♦WU'i'B'l  bl  *01 5) !♦ VI  *1  WG*Q<  5'I  *UI  bl 

♦  WII *  II I  b  I  *  1)1 6 1  I  ♦  hi:  *  I Y  I  3 »  K  i  L  I  4  HA  I  vh  X*M»*HU  I  /Al  i  I 

-4161*0141;  YIS.K, 1.7*0.  ; 

Y  1 6  ,K  i  L  I  =  I  KF*OI  3I*WA*E  16  I +00*0161  I  /  Al  31  -  I  WT*0<  4 1  ♦  A'l  6 1  *DI  6 1  »  ♦HY*IH4*lAl  11  *DI  H 
«AI2)<»I2)  )  +  WV*H5*MR*Dm  )} 

Y  1 1  *  K.  i  L  l  =  H7  *  I  V-'O*  I  )iG*C  I  2  I  +  WY*t  I  2  I-K1.4DI  1 1  I  -HIT*  A I  2  l-HF  *  Al  T I  -»  | 

Y  1 2  i  K  i  (.  I  =  H7 *  I  - V.CJ*  I  WG<- C I  I  I  •* W  Y * (' I  I  H  KL  *01  2  I  I  *HII*AI  1 1 -IIM-AT  2 1  ITGO  TO  SC? 

SOS  HS  =  V3*AI!i);  HA  =  1  IS ’•'0 1 31  ♦  AT  6)  *01 31 5  IHVAI  IT  *A  1 1 J  ♦  AI2I  *  Al  2 1 5  HE*WY*HB*H55 
Y(3iK|l.  )  =  -IA(3)*IV3*0lbl*V6*AI5I  <1)1 6  I )  ^V'AT  ;  IIF=H4’*YI  3»Ki  L  )  5 
HD*II‘>*  Yl  3i  K  i  L  I  *Y  I  3  i  K  it  I  ;  HH=WS*II0*WL*H‘)5 
Wt=VE*AI4l5  Wl  =VI*A».41  5 

Y I  4  ( K  ( L I  *  IWU*IV5*f  I+V‘/*»I4H  t«4»  »*WK*I  V5*lil  5)  *Bl  51  ♦  VITXAI  bl  *AI  5)  *01  6 1  ^01  671 

♦  YH*T  WV*V3*BI  6)  *01  M*  Vb*B  16)  *|l{  6)  *01  blXOI  ti  |  |  «V  1  *  I  Vtt+IHS*IM  5 1  ♦  AT  SJ  *1)161  ) 

♦V3*BI  b)*DI  6)  I  ♦!<(  *  I  Y  |  3  *  K,  L  I  +  HA  I  ♦riX*IIB<  HU  I  /  A  I  3  |-lilS'>  B  (4 1  ♦  Al  6  |  *U|  41  X 

Y 1 5  i  K  i  L I  *  I WL *  F, I  b  H  V.4*  Ih  B*  F 1 6  H  rtf  *11 131  )+WA*<  Vb*CI  5W9*B  I  5 1- V10*AT  51  11/ A<  3) 
♦V3*IMY*IH4*I  Al  l)*)  tl  I  4  Al  2 1 4(1(2 1  I  ♦Ht*l3 1  ?.  I  )- 1  WT*I>!  4H  Al  b )  *M  bill -Al  6X0T4IT 
YI6iK(tl‘*IWF*l)l3»  +  WH  v  |  V3*F  I  b  H  Vli*UI  b  IM'T  G  J  XhC*  I  V5*CI  6»*VV*8'1 6»  ♦  £<  6)  II /Al 
-  I  WT*l>  1 4  )  ♦HS*OI  6)  ♦  A 1 6 1  *  0 1 6  )  1  ♦  lnY*  1 1(4*  I  A'l  1  I *0  1 1 1  +  Al 2 1 *01 2 ) 0  ♦HE  *01  3 1  T  5 
Y( 1 1 K  i  L ) 

■HY* I WR*t  Vb*C 1 2 1 ♦ V9*0 1  2 1 *F 1 2 l-KK*OI 1 1 1 -MH*AI  21  -HE* A I  1 1  J'5 
YT2.K.L  )=H7*l-rfR*l  VS*t  m<V9*Dlll*tl  l  I  «KO*OT  2 1  I  *HH*  Al  H’-HF  w  Al  2)  >  5 

sc:  two;  fiio;  no  l=2  to  m-i;  on  k=i  to  n-i»  do  j=1  to  bl 

XI  J,K,U  =  X<  J,K,L»+Y(J,K*U*0IT5  fcNO;  CN05  END; 

ir  iu*2o  then  go  to  ls;  else 

IF  1 0=40  THEM  GO  TO  ls;  (  LSE 
IF  10  =  60  TllCH  GO  TO  LS;  ELSE 
IF  10=30  THEM  GO  TO  LS 5  ELSE 
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IF  10-MTT  THEN  GO  TO  LS{  ELSE  GO  TO  RIM 
LS*  DO  L*M  TO  1  BY  -II  DO  K*1  TO  NA ;  Y  1 3  ,  K  ,L  I =  X<  3,K,  D-Wl  S  Y14,K,LI* 
XI4(K,L)-VO;  DO  J  =  1 ,2,!>,6;  Yl  J,  K,  LI  «X(  J,  K  ,L)  5  ENDS  END;  END; 

DO  J-l  TO  6; 

POT  FILL  ( SYSPK I  NT )  PAGEMH)  L-M  TO  l  BY  -H  PUT  FILE  ISYSPKINT)  E 

*  D  IT  I  lDtJ(LflK«YIJ(K(LI  l!l)  K*i  TO  IfAII  I  SKI  PI2I  ,  3 1 XI4I ,  FI  3)1 .41  SKI  P  ,51  X I  *,)  ,F  I  2) 
, XI 21 ,EI 1 4,7) )  1 1 ;  END;  END; 

POT  FILE  ISYSPR1NT)  PAGE; 

Gl«0.5  DO  L*M  TO  1  BY  -I;  DO  K=1  TO  N;  ZK*KS  VI* I 2K- 1. I+DX I 
Vl-Vl/tVN-Vl 1**3; 

Yll,KfL)*X|  1,  K,LI*'XI1,K,LHXI2,K,L»*XI2,K,L»;  Yl  2  ( K,  L )  *V1  ♦  Yl  l ,  K(  L 1 1  END; 

DO  K* 1  TO  N-2  BY  25 

G1IL)=G1IL)*YI2,K,LM',.*YI2,K+1,L)*YC2,K*2,L){  END; 

Gl( L  )*3< 141 593E1C/WG*G1 (L )*DX{ 

END!  DO  L*M  TO  1  BY  -I;  PWT*0. 1  00  LL*L  TO  H-2  BY  25  PW1  -I’WT  <G1 1 LL ) 

♦  4.*G1(LL+1)«-G1(LL*2)S  ENDS  I F  LL<K-2  THEN  PWI =PWT  + . 5* I  GllH-ll+Gl IM) I ! 
PWT*DY*PWT/3. ;  DO  K*1  TO  NAS  Yl  1 ,K, l )*SQRTI Yl 1 ,K ,L 1 1  5 END; 

PUT  FILE  I SYSPRI NT I  L 

DIT  I  lDfJiLf  IKiYIlrKiLI  DO  K=1  TONAI)  I SK IPC  21 , 3 1 XI 41  (FI  3 )  I  (41  SK I P  ,  51  X I !»l ,  F I  2) 
(XI 2). £114,71111; 

PUT  FILE  ISYSPKINT)  EDIT  IG1 ( L ) , PhT I  ISK I P, 21 X I  21 , E I  14 , 7 ))) S  ENDS 
KK:  end;  IF  ic=l  Then  do;  ic*o ;  ep*3.E3;  IJ-100;  ZIJ-IJ;  CHG*1.; 

D1T-1.E-7; 

GO  TO  RT;  end;  end; 
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ILLUSTRATIONS 


Figure  1.  Fore  and  Aft  Symmetry  of  the  Laser  Pulse. 

The  on -axis  magnitude  of  the  electric  field,  |e|,  is  shown  as  a 
function  of  |z  -  zc  J  at  four  different  times.  As  shown  in  Eq.  (200), 
zc  w  ct  locates  the  center  of  the  pulse.  The  solid  curves  show  the  lead¬ 
ing  edge  of  the  pulse  and  the  dotted  curves  depict  the  trailing  edge.  Dis¬ 
tances  along  the  z-axis  are  expressed  in  units  of  the  grid  size:  Az  =  0.45km. 
The  quantity  JeJ  shown  is  defined  in  Eq.  (197),  so  that  the  exponential 
damping  factor  is  not  included  in  the  graphs. 


Figure  2.  Details  of  the  z-profile  of  the  Electric  Field. 

For  detailed  comparison,  the  z-profile  of  the  on-axis  electric  field,  |e|, 

is  shown  at  several  times.  The  unit  Az  is  used  for  distances  along  the 

z-axis.  The  curves  have  been  displaced  to  the  left  and  the  leading  edges 

made  to  coincide  at  height  1000  for  |e|.  The  abscissa  for  this 

\cmJ/ 


intersection  of  the  curves  has  been  labeled  42.  5,  the  location  of  this  point 
at  t  =  0  . 


Figure  3.  Radial  Profiles  of  the  Laser  Pulse. 

J  j  _5 

The  radial  profile  of  |E|  is  shown  for  t  =  0  and  for  t  =  10  sec  for 

slices  taken  through  the  on-axis  maximum,  Zp,  in  the  z-profile. 

Figure  4.  Details  of  the  Radial  Profile  of  the  Laser  Pulse. 

Details  of  the  radial  profile  are  shown  for  various  times.  In  all  four  cases 
the  radial  slice  through  the  on-axis  maximum  of  the  z-profile  is  exhibited. 
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Figure  5.  Off-axis  Maximum  in  Laser  Pulse. 

An  off-axis  maximum  in  the  radial  profile  is  shown  at  t  =  10  ^sec.  The 
slice  shown  exhibits  the  radial  profile  at  z  =  31,  whereas  the  principal 
peak  of  the  pulse  is  on-axis  at  z  =  35.  3  .  The  slice  at  z  =  31  contains 
the  greatest  off-axis  effect  and,  therefore,  locates  the  two  secondary  peaks 
which  have  developed  in  the  pulse.  These  secondary  peaks  are  also  indi¬ 
cated  on  Figure  9. 


Figure  6.  Phase  of  the  Electric  Field  Amplitude. 

Phase  information  is  presented  by  showing  |E.  |  as  a  function  of  z  at 
/  * 

t  =  6  x  10”  sec.  For  comparison,  the  dashed  curve  shows  |e|.  Equa¬ 
tions  ( 162  )  and  (197  )  of  the  text  define  JEjJ  and  |e|. 

Figure  7.  Phase  of  the  Electric  Field  Amplitude. 

At  t  =  8  X  10 "  sec. ,  |EjJ  and  |e|  are  shown  on -axis  as  functions  of  z. 
Two  nodes  have  developed  and  E^  is  negative  in  the  region  of  the  power 
peak.  The  sign  of  E^  in  the  various  regions  is  indicated  on  the  figure. 
The  node 8  are  also  shown  in  Figure  9- 


Figure  8.  Phase  of  the  Electric  Field  Amplitude. 

At  t  =  10'^sec.  ,  )EjI  and  |e|  are  shown  on-axis  as  functions  of  z. 
The  sign  of  Ej  is  indicated  in  the  various  regions.  There  are  now  four 
nodes.  The  nodal  curves  are  phase  fronts  and  are  shown  in  detail  in 
Figure  9. 


Figure  9.  Configuration  of  the  Laser  Pulse. 

Various  properties  of  the  pulse  are  shown  in  the  rz -plane.  The  loca¬ 
tion  of  the  peak  in  the  z -profile  is  shown  as  a  function  of  r  at  t  =  0  and  at 

-5  -5 

t  =  10  sec.  The  phase  fronts  with  E,  =  0  are  shown  at  t  =  10  sec.  The 

*  -6 

open  circles  locate  the  z -profile  nodes  of  E,  at  t  =  8  X  10  sec.  The 

1  _5 

small  squares  locate  the  secondary  maxima  of  the  pulse  at  t  =  10  sec. 
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Figure  10. 


Figure  11. 


Figure  12. 


Fig  .re  13. 


Figure  14. 


Figure  15. 


Figure  16. 


On-axia  Temperature  Distribution. 

The  on-axis  temperature  increment,  (T  -  TQ),  is  shown  as  a  function  of 
z  at  t  =  10'^sec. 

Radial  Temperature  Profile. 

-5 

The  radial  profile  of  the  temperature  increment  is  shown  at  t  =  10  sec. 
for  the  slice  through  the  maximum  of  the  z -profile.  This  maximum  is  at 
z  =  33  as  may  be  seen  in  Figure  10. 

Fluid  Velocity  Distribution, On-axis. 

The  z -component,  v  ,  of  the  fluid  velocity  is  shown  on-axis  as  a  function 
-5  z 

of  z  at  t  =  10  sec.  A  double  log  plot  is  used  which  omits  values  of  v 

-3  -3  z 

between  10  cm/sec  and  -10  cm/sec. 

Radial  Fluid  Velocity  Distribution. 

The  radial  component,  v  ,  of  the  fluid  velocity  is  shown  as  a  function  of 
-5  r 

r  at  t  =  10  sec.  The  slice  is  taken  at  z  =  32,  the  location  of  the  "center 
of  velocity"  shown  in  Figure  12. 

Fluid  Density  Distribution,  On-axis. 

The  on-axis  fluid  density  decrement,  -(p-p  ),  is  shown  as  a  function  of 

-5  ° 

z  at  t  =  10  sec. 

Radial  Fluid  Density  Distribution. 

-5 

The  radial  density  distribution  is  exhibited  as  a  function  of  r  at  t  =  10  sec. 
The  slice  is  taken  at  z  =  32,  the  location  of  the  density  minimum  detailed 
in  Figure  14.  A  double  log  plot  is  used  which  omits  values  between 

m-’-enj  lnd  . 

cm  cm 

Parade  of  Laser  and  Fluid  Pulses,  On-axis. 

On  a  double  log  plot,  the  various  laser  and  fluid  variables  are  simultaneously 

-5 

plotted  versus  z  at  t  =  10  sec  so  that  the  spatial  location  of  the  various 
pulses  can  be  visualized. 
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Figure  17.  Locations  and  Widths  of  Laser  and  Fluid  Pulses. 

The  location  and  full  — -widths  of  the  various  laser  and  fluid  pulses  are 

e 

shown  versus  z.  The  peak  to  valley  distance  is  shown  for  v^.  The  initial 
and  final  locations  of  the  laser  power  peak  are  also  indicated. 

Figure  18.  Self-curving  Effect  of  a  Laser  Beam  in  the  Presence  of  a  Wird . 

The  self-curving  effect  of  a  laser  beam  in  the  presence  of  a  wind  is  in  - 
dicated.  The  angle  the  deflection  angle,  increases  with  z.  The  density, 
permittivity,  or  index  of  refraction  profile  is  sketched  to  indicate  the  cause 
of  the  curving. 


Figure  19.  Thermal  Blooming  of  a  Laser  Beam. 

The  thermal  blooming  of  a  laser  beam  in  the  absence  of  wind  is  indicated. 

The  angle  the  deflection  angle  of  the  outer  edge  of  the  beam,  increases 
with  z.  The  density,  permittivity,  or  index  of  refraction  profile  is 
sketched  to  indicate  the  cause  of  the  blooming. 

Figure  20.  Cross -section  for  Scattering  into  Various  Curvatures. 

The  normalized  cross-section  for  scattering  of  power  into  various  curva¬ 
tures  is  sketched  for  small  wind  speeds,  X  =  0,1,3.  The  distribution  ex¬ 
tends  from  to  *or  each  X.  Already  at  X  =  3,  the  dis¬ 
tribution  is  sharply  peaked  at  ?  =  ?  .  All  the  curves  must  pass 

tn&x 

through  the  vertical  axis  at  height  0. 5  and  have  an  asymptote  at  5/St,  =  2  . 

D 

Figure  21.  Effective  Optical  Elements  for  Blooming  . 

The  effects  of  a  wind  on  the  beam  are  indicated.  When  X  =  0,  pure  ther¬ 
mal  blooming  occurs;  the  light  rays  bend  as  if  passing  through  a  symmetrical 
double  concave  lens  at  each  point  in  the  path.  For  X  t  0  the  effective  lens 
is  converted  into  a  large  prism  with  a  small  defocusing  lens  perched  on 
the  downwind  vertex  of  the  prism. 
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Figure  22.  Peak  Displacement  at  Detector  for  Experiment  of  Smith  and  Gebhardt. 

For  the  experiment  described  in  the  text,  the  beam  is  viewed  with  a  detec¬ 
tor  located  20cm  from  the  end  of  a  100cm  wind  tunnel.  The  peak  displace¬ 
ment  into  the  wind,  Y^/a,  is  shown  as  a  function  of  velocity.  The  nota¬ 
tion  is  explained  in  Part  VI  of  the  text. 


Figure  23.  Intensity  Profiles  at  Detector  for  Experiment  of  Smith  and  Gebhardt. 

Intensity  profiles  predicted  by  equation  (290)  of  the  text  are  shown  at 
aL  =  0.5  for  vq  =  2cm/sec  and  lOcm/sec.  The  original  beam  is  also 
shown. 
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(lb)  t  =  6x10  sec  (lc)  t  =  8x10  sec  I  (Id)  t  =  10  sec 


Figure  1.  Fore  and  Aft  Symmetry  of  the  Laser  Pulse 
Leading  edge:  solid  curves. 

Trailing  edge:  dashed  curves 


10  20  30  40  |  50 

42.  5 


(A  z  =  0.  45  Ion) _ > 

Figure  2.  Details  of  the  z-profile  of  the  Electric  Field 
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Figure  11.  Radial  Temperature  Profile 
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Figure  12  Fluid  Velocity  Distribution  On- Axis 
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Figure  13.  Radial  Fluid  Velocity  Distribution 
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Figure  17 


Locations  and  Widths  of  Laser  and  Fluid  Pulses 


Figure  18.  Self-curving  Effect  of  a  Laser  Beam  in  the  Presence  of  a  Wind 
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Figure  21.  Effective  Optical  Elements  for  Blooming 
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Figure  23.  Intensity  Profiles  at  Detector  for  Experiment  of  Smith  and  Gebhardt 
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19.  ABSTRACT 


In  this  work  instabilities  in  tbo  system  of  equations  describing  electromagnetic 
wave  propagation  and  fluid  dynamics  are  discussed.  The  basic  equations  are 
discussed  in  Parti,  and  in-Part-H  the  linearized  stability  analysis  is  presented 
along  with  an  evaluation  of  the  threshold  for  these  growing  waves.  To  follow  the 
growth  of  the  disturbances,  computer  studies  were  undertaken.  In  the  course  of 
these  studies  it  became  apparent  that  there  was  some  merit  to  introducing  a  new 
concept  to  judge  the  value  of  an  algorithm  for  computing  the  solutions  of  a  system 
of  partial  differential  equations.  This  concept  was  called  "utility",  and  is  discussed 
in-Part-UI*  along  with  several  examples  of  its  application  to  simpler  partial  differ¬ 
ential  equations^  The  advantage  of  this  concept  is  that  it  is  relatively  easy  to  apply 
to  complicated  systems  of  partial,. differential  equations,  whereas  the  stability 
concept  leads  to  a  very  complicated  procedurq  for  deciding  on  the  value  of  a 
numerical  routine.  In  Part  IV  are^resented-uie  results  of  a  calculation  of  beam 
distortion  for  a  very  high  intensity  pulse  propagating  through  air  for  several 
kilometers.  Analytical  arguments  are  advanced  which  suggest  that  tb*  qualitative 
features  of  the  distortions  are  correct,  which  lends  credence  to  the  computer  output, 
Speed  and  memory  size  in  a  computer  place  certain  restrictions  on  one's  ability 
to  investigate  phenomena  in  the  laser  beam  problenSr>  In  the  attempt  to  calculate 
distortions  of  the  type  predicted  by  the  linearized  instability  analysis,  cylindrical 
symmetry  was  imposed  on  the  problem  in  order  to  faciUxte  the  computer  calculatioi 
Had  this  not  been  necessary,  or  had  some  other  independent  variable  been  eliminate^ 
(continued  on  next  page) 
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13,  Abstract  (continued) 

rather  than  the  angle  about  the  beam  axis,  much  more  pronounced  evidence  of 
beam  and  fluid  instabilities  would  likely  have  been  observed  for  substantially 
lower  powers,  powers  that  maybe  achievable.  Arguments  supporting  this 
proposition  are  contained  in  Part  V. 

Part  VI  of'thjjs  report  contains  an  analytic  discussion  of  beam  bending  and 
the'rmal  blooming  for  a  slab  beam  propagating  through  a  wind.  A  formula  is 
derived  which  provides  for  the  transition  between  two  regimes  in  which  con¬ 
duction  and  forced  convection,  respectively,  dominate  the  dissipation  of  heat 


deposited  in  the  medium  from  the  laser  beam., 
useful  for  the  analysis  of  several  experiments. 


This  formula  appears  to  be 


